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I.    INTRODUCTION 

A.  PROBLEM  DEFINITION 

In  recent  years  considerable  effort  has  been  put  into  the 
improvement  of  infrared  sensors.  Increased  sensor  resolution 
was  one  of  the  more  ardently  pursued  goals.  The  drive  behind 
the  desire  for  better  resolution  came  from  several  areas.  One 
of  the  driving  forces  was  the  military's  need  for  infrared 
sensors  that  were  capable  of  detecting  heat  signatured  objects 
at  long  range. 

The  specific  concern  in  this  thesis  is  the  infrared 
detection  of  heat-signatured  objects  with  sensors  at  altitudes 
of  several  hundred  kilometers.  Specifically,  the  thesis  will 
be  concerned  with  the  detection  and  tracking  of  airborne 
targets  such  as  medium  range  missiles  and  aircraft.  The 
detection  and  tracking  problem  will  be  addressed  from  a 
generalized  infrared  sensor  viewpoint.  This  is  done  with  the 
desire  to  produce  methods  with  the  widest  applicability. 

B.  BACKGROUND 

The  signature  size  and  velocity  of  the  targets  have 
considerable  bearing  on  the  direction  taken  in  solving  this 
problem.  Typical  sensor  geometric  resolutions  are  on  the 
order  of  .08  mrad.   [Refs.  1  through  5]    This  produces  a 


footprint  on  the  earth's  surface  of  approximately  48  m  for  a 
sensor  at  600  km  (1  /irad=l  m  at  1000  km)  .  The  type  of  targets 
considered  in  development  of  the  problem  have  signatures 
smaller  than  the  footprint  of  the  sensor.  This  will  be 
assumed  to  hold  true  even  with  the  effects  of  point  spreading 
and  optical  blur  inherent  in  the  sensors  [Refs.  3  and  4]. 
Sensor  sampling  rates  are  typically  on  the  order  of  1  sample 
per  resolution  cell,  here  .08  mrad  [Ref.  3].  The  resultant 
intensity  of  a  sample  containing  a  target  is  due  to  the  sum 
of  target  thermal  intensity  and  the  background  intensity 
within  the  footprint  [Ref.  2].  Thus  the  targets,  with  a 
modified  intensity,  in  a  discrete  image  are  nominally  about 
one  pixel  in  dimension. 

The  velocity  of  the  targets  is  considerably  larger  than 
the  drift  rate  due  to  sensor  motion  (e.g. ,  two  percent  of  the 
detector  footprint)  and  jitter  [Ref.  4].  If  the  image  scan 
rate  is  typically  30  image  frames  per  second,  normal  target 
movement  between  images  will  be  less  than  two  pixels  per 
frame.  A  majority  of  the  targets  considered  would  actually 
move  at  less  than  one  pixel  per  frame. 

Many  of  the  infrared  systems  contain  more  than  one  focal 
plane  assembly  operating  at  different  wavelengths.  The 
resulting  images  can  be  dealt  with  as  individual  images  each 
with  its  own  statistical  characterization.  In  a  single  focal 
plane  the  scenes  may  be  characterized  by  several  2-D  random 


processes.  The  transition  from  one  process  to  an  adjacent 
process  may  be  correlated  or  uncorrelated.  The  processes 
describe  such  items  as  terrain,  and  multiple  layers  of  clouds 
and  have  a  mean  value  dependent  on  the  scene  and  the  spectral 
band.  [Ref.  4]  Additionally  the  effects  due  to  atmospheric 
jitter,  background  scintillation,  and  electronic  noise  may  be 
modeled  by  a  combination  of  white  and  correlated  gaussian 
noise. 

C.   CONCEPT  AND  CONTENT  OF  THE  THESIS 

The  preceding  sections  have  provided  the  background  to 
develop  a  conceptual  solution  for  the  problem.  It  can  be  seen 
that  ultimately  the  detection  algorithm  will  need  to  detect 
pixel-size  targets  with  intensities  that  may  differ  only 
slightly  from  the  surroundings.  In  tracking,  the  primary 
region  of  operation  is  seen  to  involve  target  motion  of  one 
to  two  pixels  between  successive  image  frames. 

In  this  thesis  the  solutions  will  be  developed  for  images 
in  a  single  focal  plane.  Their  applicability  will  be  general 
in  nature  and  do  not  depend  on  the  frequency  band  or  the 
detailed  sensor  characteristics.  Some  ideas  for  merging  the 
results  from  applying  the  algorithms  in  several  focal  planes 
simultaneously  will  be  addressed  in  the  final  chapter. 

The  next  four  chapters  deal  with  the  detection  of  the 
target  and  Chapter  VI  deals  with  the  tracking.   Initially  a 


method  will  be  proposed  for  modelling  the  background  processes 
present  in  the  image  (terrain,  clouds,  etc.).  In  subsequent 
developments  the  background  processes  will  be  assumed  to  be 
homogeneous.  Surveys  of  the  literature  show  that  frequently 
90%  of  the  image  may  be  characterized  by  a  single  random 
process.  The  modelling  process  will  be  performed  with  fixed 
length  spatially  adaptive  two-dimensional  prediction  filters. 
Adaptive  filters  used  in  the  modelling  process  are  the  topics 
of  Chapters  III  and  IV.  In  these  chapters  the  two-dimensional 
Least  Mean  Square  (LMS)  and  Recursive  Least  Squares  (RLS) 
filters  are  derived  and  tested. 

The  detection  process  is  the  subject  of  Chapter  V.  Linear 
prediction  produces  an  error  residual  process  which  is  stored. 
This  error  process  contains  the  unpredictable  information  in 
the  image.  In  Chapter  V  a  method  based  on  statistical 
significance  testing  will  be  used  to  separate  the  error 
residual  into  levels  of  predictability.  The  error  residual 
present  at  the  target  location  will  be  shown  to  be  anomalous 
by  comparison.  This  results  in  highlighting  the  targets  after 
thresholding  to  achieve  a  low  constant  false  alarm  rate. 

In  Chapter  VI  a  point  target  tracker  is  developed  and 
tested.  The  tracker  is  the  based  on  a  maximum  a  posteriori 
(MAP)  estimate  of  the  target  position.  In  the  overall  system 
the  detection  algorithm  passes  its  results  to  an  autonomous 
tracking  algorithm.   Consideration  of  the  logical  interface 


between  the  detection  algorithm  and  the  tracking  algorithm  are 
addressed  in  Chapter  VI.  The  precept  throughout  is  the 
development  of  a  set  of  algorithms  that  is  reasonably  generic. 
Algorithms  developed  extend  readily  to  images  other  than 
infrared. 


II.  MODEL 

A.   GENERAL  DESCRIPTION 

In  order  to  pursue  the  objectives  presented  in  the 

previous  chapter  it  is  necessary  to  establish  a  framework  from 
which  to  work.  Throughout  this  thesis  the  image  to  be 
processed  is  taken  to  be  comprised  of  three  components.  The 
three  components  maybe  viewed  as  a  background  generation 
process,  a  corrupting  noise  process,  and  a  target  generation 
process.  Since  target  detection  is  the  ultimate  objective, 
it  will  be  necessary  to  separate  the  three  components. 

Background  generation  processes  have  been  studied  by  many 
researchers;  the  result  is  a  wide  variety  of  models.  If  the 
image  can  be  divided  into  statistically  similar  regions,  then 
each  region  may  be  modeled  on  an  individual  basis.  It  will 
further  be  assumed  that  although  the  overall  image  may  be 
nonstationary,  the  statistics  in  the  separate  regions  are 
individually  stationary.  The  region  within  the  image  is  taken 
to  be  the  result  of  a  linear  shift  invariant  operation  on  a 

wide  sense  stationary  Gaussian  distributed  white  noise  source 

2 
of  variance  o   .   [Ref.  6]   A  block  diagram  of  the  2-D  linear 

operation  is  shown  in  Figure  2-1. 


Figure  2-1.   Block  Diagram 


The   particular   linear   operation   to   be   used   is   a   two 

dimensional  autoregression  operation  which  may  be  written  as: 

x(n,m)  =  E  E  h(£,k)  x(n-£,m-k)  +  w(n,m) 
I   k 
(«,k)  e  (3 

(2-1) 

where  /J  is  some  chosen  region  of  support  for  the  filter 

coefficients.   This  model  has  been  termed  the  White  Noise 

Driven  Representation  (WNDR)  and  the  transfer  function  may  be 

viewed  as  an  all-pole  IIR  filter  with  a  constant  numerator 

driven  by  a  white  noise  source.   The  filter  is  said  to  be 

"causual"  (in  the  2-D  sense)  if  the  computation  of  x(n,m)  does 

not  require  knowledge  of  "future"  values  of  the  process. 

Otherwise  it  is  said  to  be  noncausal.   The  definition  of 

causality  is  tied  to  the  scan  pattern  used  in  processing  the 

image.   Figure  2-2  depicts  the  issue  of  causality  for  a  row 

scan  direction. 


-^  scan    direction 


Figure  2-2.  Casuality  Depiction 


If  it  is  possible  to  determine  the  coefficients  which  describe 
the  regional  process,  it  then  becomes  possible  to  inverse 
filter  the  region.  This  will  ideally  leave  only  the  driving 
noise  process. 

Noise  present  in  the  image  is  assumed  to  be  additive 
although  not  necessarily  gaussian  or  white.  If  the  noise  is 
colored,  the  same  procedure  used  to  identify  and  eliminate  the 
regional  process  would  be  utilized  on  the  noise.  The  result, 
as  before,  is  a  white  noise  process. 

The  most  difficult  of  the  three  components  of  the  image 
to  characterize  is  the  target  process.  If  a  specific  model 
for  the  target  is  specified  then  the  possibility  exists  that 
the  overall  system  would  be  too  restrictive.   In  order  to 


avoid  this  pitfall  it  is  proposed  that  the  only  identifying 
features  of  the  target  are  its  pixel/sub-pixel  size, 
opaqueness,  and  non-zero  mean. 

B.   PREDICTION  ERROR  FILTER 

It  has  been  shown  [Refs.  6  and  7],  if  the  regional  process 

is  stationary  and  under  specific  conditions  of  process 
support,  that  an  optimal  inverse  filter  solution,  in  the  MMSE 
sense,  for  the  original  process  may  be  found.  The  drawbacks 
to  the  proposed  solution  are  that  the  filter  support  with  few 
exceptions  will  be  of  infinite  extent.  Further  a  priori 
knowledge  of  the  2-D  correlation  function  or  power  spectrum 
and  possibly  the  use  of  spectral  factorization  techniques  are 
required  to  find  the  filter.  Physical  realization  of  the 
filter  would,  in  any  event,  dictate  truncation  of  the 
resultant  filter  support  region  prior  to  implementation.  Once 
the  filter  support  region  has  been  truncated,  it  no  longer  is 
capable  of  matching  the  original  regional  autocorrelation 
function  and  may  not  be  minimum  phase  and  therefore  not 
stable.   [Refs.  8  and  9] 

The  above  considerations  prompted  the  use  of  reduced  order 
models  for  the  regional  processes  and  adaptive  techniques  for 
fitting  the  reduced  order  model.  The  support  of  the  filter 
coefficients  of  the  reduced  order  model  is  causal  and  has 
finite  support  as  shown  in  Figure  2-3. 
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Figure  2-3.   Reduced  Order  Model  Filter 
Support 


The  resultant  reduced  order  model  is 

x(n,m)  =  S  S  h(£,k)  x(n-£,m-k)  +  v(n,m) 
H   k 
(£,k)  e  a 

(2-2) 

where  v(m,n)  is  a  noise  source,  not  necessarily  white,  that 

would  be  required  to  exactly  reproduce  the  original  regional 

process. 

The  driving  noise  process  may  also  be  viewed  as  the  error 

residual  in  the  process  produced  by  the  inverse  of  the  reduced 

order  model  operating  on  the  original  process  as  shown  below. 
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e(n,m)  =  x(n,m)  -  E  Z  h(£,k)  x(n-£,m-k) 

I   k 
(£,k)  e  a 

(2-3) 

This  representation  is  the  prediction  error  filter  (PEF)  form 

[Ref.   9].     In  practice   if  the   support  of  the   filter 

coefficients  in  the  previous  model  is  causal,  and  sufficiently 

large,   the  process  complex  spectral  density  function  is 

analytic  in  the  neighborhood  of  the  unit  bicircle,  and  a 

minimum  error  variance  fitting  technique  is  employed  to  fit 

the  PEF  to  the  process;  then  the  spectral  density  function 

can  be  closely  matched.   In  this  situation,  the  noise  term 

e(m,n)  is  very  nearly  white  and  gaussian  if  the  original 

process  is  gaussian. 

In  practice  the  image  background  can  only  be  assumed 

homogeneous    (i.e.,    stationary)    in   a   small   region. 

Consequently  the  model  will  have  to  be  fitted  to  each  local 

region  or  adapted  during  the  image  processing.   Figure  2-4 

depicts  a  section  of  an  image  used  for  this  processing.   The 

point  x(n,m)  is  the  point  at  which  the  model  of  the  process 

as  presented  in  (2-3)  is  to  be  fitted.   The  prediction  window 

represents  the  nonzero  support  of  the  summation  term  in  (2-3) 

which  will  be  referred  to  as  the  predictor.   The  prediction 

error  filter  is  defined  with  a  finite  rectangular 
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Figure  2-4.   Reduced  Order  Model  Fitting 


support  of  size  PxQ  and  for  convenience  of  notation  the 
coefficients  are  assigned  as  follows 


a(0,0)  =  1 


(«,k)  =  (ofo) 


a(£,k)  =  -h(£,k)   0  <  i   <   P-l, 

0  <  k  <  Q-l 

(2-4) 

The  points  underlying  the  prediction  error  filter  mask  may  be 

concatenated  into  a  column  vector 
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X  = 


x(n,m) 
x(n-l,m) 


x(n-P+l,m) 

x(n,m-l) 
x(n-l,m-l) 


x(n,m-Q+l) 
x(n-l,m-Q+l) 


x(n-P+l,m-Q+l) 

(2-5) 
The  point  ordering  establishes  the  past  (i.e.,  causality)  with 
respect  to  the  point  x(n,m) .  From  this  vector  the  correlation 
matrix  can  be  defined  as 

R  =  E{x  xT} 

(2-6) 
Since  the  region  is  assumed  stationary,  the  correlation  matrix 
is  block  Toeplitz.   The  optimal  prediction  coefficients  are 
determined  by   solving  the  Normal   Equations   for   linear 

prediction  as  shown  below: 

2  . 
R  a  =  ai 

x —     e  — 

(2-7) 

where  i  is  a  unit  vector  that  is  equal  to  one  in  the  first  row 

2  .  . 

and  zero  elsewhere.   The  term  a     represents  the  minimum  error 

variance,  in  the  MMSE  sense,  of  the  white  residual  error  of 
the  PEF  using  the  optimal  coefficients.   In  this  thesis  the 
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2-D  Least  Mean  Square  (LMS)  algorithm  is  used  to  obtain  an 

approximate  solution  to  (2-7)  in  a  local  area. 

Under  the  assumption  of  local  stationarity  it  is  also 

possible   to   produce   a   set   of   prediction   error   filter 

coefficients  for  the  region  directly  from  the  data.    The 

method  employed  to  produce  that  solution  is  to  construct  a 

family  of  equations  of  the  form  of  (2-3)  which  represent  the 

realization  of  the  prediction  error  filter  at  each  point  in 

the  region.   These  equations  are  then  solved  for  the  filter 

coefficients  subject  to  the  constraint  that  the  sum  of  their 

error  variances  is  minimized.   The  quantity  to  be  minimized 

is: 

J(x,a)  =  E{2  S  [x(n,m)  -  x(n,m)]^} 
n  m 
(n,m)  e  R 

(2-8) 

where  x(n,m)  is  the  prediction  estimate  (2-9)  for  the  point 

at  position  (n,m)  in  region  R,  which  has  been  previously 

defined  as  a  region  in  which  the  process  may  be  considered 

stationary,  given  by 

x(n,m)  =  2  S  h(£,k)  x(n-£,m-k) 
£  k 

(£,k)  e  a 
(£,k)  f    (0,0) 

(2-9) 

This  procedure,  known  as  "least  squares,"  has  been  employed 

in  many  applications.   In  particular  Ref.  10  has  specific 
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bearing  on  this  thesis.  The  recursive  form  of  the  least 
squares  (RLS)  technique  has  been  used  in  the  thesis. 

The  final  result  of  either  the  LMS  or  the  RLS  method  is 
a  "whitening"  of  the  background  process.  The  details  of  the 
actual  processing  will  be  covered  in  the  remainder  of  the 
thesis. 

C.   SIGNIFICANCE  TESTING 

Under  hypothesis  testing  a  threshold  is  established  which 

delineates  the  boundary  between  two  hypothesises  Hq  and  H,. 
This  can  be  seen  to  be  a  two  sided  testing  procedure.  The 
threshold  is  established  by  choosing  the  accceptable 
probability  of  rejecting  Hq  when  Hq  is  true,  P(type  I  error) 
and  minimizing  the  probability  of  choosing  Hq  when  H-.  is  true, 
P(type  II  error) .  This  method  presupposes  known  distributions 
for  the  two  hypotheses.  It  has  been  stated  previously  that 
a  robust  target  detection  system  should  minimize  the 
specifications  on  the  target  process.  Hypothesis  testing 
requires  too  much  to  be  known  about  the  target. 

Significance  testing  is  a  one  sided  testing  process  with 
only  one  hypothesis  Hq,  refered  to  as  the  "null  hypothesis." 
The  procedure  is  to  completely  specify  the  probability 
distribution  of  Hq,  then  compute  a  test  statistic  from  the 
event  to  be  tested  against  Hq.  The  test  statistic  is  a 
function  of  the  event  and  is  based  solely  on  the  event.   The 
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decision  to  accept  or  reject  Hq  after  the  statistic  is 
compared  to  Hq  is  subjective,  in  that  the  establishment  of 
boundary  between  the  two  decisions  is  at  the  descretion  of  the 
designer.  Using  the  test  statistic  and  the  hypothesized 
probability  distribution  of  Hq  a  region  of  values  for  the  test 
statistic  is  determined  where  Hq  will  be  rejected.  This 
determines  the  probability  of  a  type  I  error,  which  is  known 
as  the  significance  level  of  the  test.  The  significance 
level  may  be  based  on  desires  to  minimize  the  false  alarm  rate 
(time  between  type  I  errors) ,  maintain  a  constant  false  alarm 
rate,  etc.   [Ref.  6] 

The  significance  testing  used  for  target  detection  in  this 
thesis  involves  a  hypothesized  gaussian  distributed  white 
noise  background  (Hq) .  This  noise  is  ideally  the  result  of 
decorrelating  the  background  as  previously  mentioned.  The 
specifics  of  the  testing  will  be  covered  in  Chapter  V. 
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III.     LMS  FOR  TOO-DIMENSIONAL  PREDICTION 

A.   THEORY 

In  the  previous  chapter  the  use  of  a  PEF  with  a  reduced 
order  predictor  model  was  proposed  as  a  method  to  decorrelate 
the  background  process.  This,  however,  was  contingent  on 
choosing  the  predictor  coefficients  such  that  the  error 
variance  at  the  output  of  the  PEF  was  minimized,  see  Figure 
3-1. 
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Figure  3-1.   Prediction  Error  Filter 
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This  chapter  extends  the  traditional  LMS  algorithm  to  the 

problem  of  two-dimensional  PEF.  In  this  form  the  coefficients 

of  the  predictor  are  adaptively  changed  to  minimize  the  mean 

squared  error.    [Refs.   11  and  12]    The  criterion  to  be 

minimized  is 

£  =  E[e(n,m;i)2] 

(3-1) 

The  error  term  e(n,m;i),  for  the  predictor  mask  of  size  PxQ 

with  first  quadrant  support  shown  in  Figure  3-1,  is  defined 

as 

P-l  Q-l 
e(n,m;i)  =  x(n,m)  -   2   2  h(£,k;i)  x(n-£,m-k) 

2=0  k=0 
(«,k)  t    (0,0) 

(3-2) 

where  the  index  i  represents  the  ith  update  of  the  filter.   If 

the  error  is  squared  the  criterion  becomes 

P-l  Q-l 
£   =  E{x(n,m)2}  -  2E{x(n,m)   S   S  h(£,k;i)  x(n-£,m-k)} 

1=0   k=0 
P-l  Q-l  P-l  Q-l 
+E{  I   S   S   2  h(g,h;i)  h(£,k;i)  x(n-g,m-h)  x(n-£,m-k)} 
g=0  h=0  ^=0  k=0 
(g,h)  f    (0,0) 
(£,k)  f    (0,0) 

(3-3) 

The  signal  x(n,m)  is  assumed  to  be  stationary  and  ergodic. 

It  will  be  shown  later  that  the  stationarity  assumption  can 

be  relaxed  and  still  produce  adequate  results.   The  result  of 
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(3-3)  is  a  quadratic  error  surface  that  has  a  unique  minimum. 

Differentiating   (3-3)   with   respect   to   the   predictor 

coefficients  yields 

P-l  Q-l 

d£ =  v(i)  =  -2E{[x(n,m)  -  2   2  h(g,h;i)  x(n-g,m-h)]- 

ah(£,k;i)  g=0  h=0 

x(n-£,m-k) ) 
o  <  £  <  p-i 
0  <  k  <  Q-l 
(£,k)  t    (0,0) 

(3-4) 

This  is  the  gradient  of  the  error  surface  at  update  i  mapped 

into  the  predictor  coefficient  space.   The  term  x(n-£,m-k) 

represents  the  array  of  image  points  under  the  prediction 

mask;  the  gradient  v(i)  also  is  in  array  form.  At  the  minimum 

point  on  the  surface  the  gradient  vanishes.   The  result  of 

minimizing  (3-4)  and  substituting  into  (3-3)  produces  the 

following  expression  for  the  minimum  mean  square  error  (MMSE) : 

P-l  Q-l 
Cmin  =  E{x(n,m)^}  -  E{x(n,m)   2   2  h(£,k;i)  x(n-£,m-k) } 

£=0  k=0 

(3-5) 

B.   THE  STEEPEST  DESCENT  ALGORITHM  (SDA) 

Central  to  the  theme  of  the  LMS  algorithm  is  the  method 

of  steepest  descent.  This  exploits  the  quadratic  character 
of  the  error  surface  by  iterative  adjustment  of  the  filter 
coefficients.   The  coefficients  are  adjusted  by  an  amount 
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proportional  and  opposite  in  direction  to  the  gradient  at  the 

current  location  on  the  error  surface,  that  is 

h(£,k;i+l)  =  h(£,k;i)  +  /i(-V(i)) 

0  <  £  <  P-l 
0  <  k  <  Q-l 
(!,k)  t    (0,0) 

(3-6) 

Here,  as  in  the  one  dimensional  LMS  [Ref.  11],  the  gradient 

is  approximated  by  removing  the  expectation  operator  from  (3- 

4)  and  recognizing  the  bracketed  term  as  the  prediction  error, 

e(n,m;i),  at  update  i.  The  form  of  the  coefficient  update  (3- 

6)  now  becomes 

h(£,k;i+l)  =  h(£,k;i)  +  2/x{e(n,m;i)  x(n-£,m-k)} 
0  <  i    <    P-l 
0  <  k  <  Q-l 
(£,k)  t     (0,0) 

(3-7) 

The  term  /x  is  a  scalar  known  as  the  step  size  and  controls  the 

rate  of  adaption  [Ref.  11] . 

C.   ALGORITHM  ANALYSIS 

In  order  to  implement  the  LMS  algorithm  it  necessary  to 

use  (3-2)  and  (3-7) ,  and  to  choose  an  appropriate  value  for 
/i  and  an  initial  condition  on  the  coefficient  array. 

The  LMS  algorithm  then  attempts  to  track  the  error 
surface,  under  the  assumption  of  stationarity,  by  adjusting 
the  filter  coefficients.   The  result  is  an  error  surface  that 
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is   inherently   changing   even   if   the   original   image   is 
stationary. 

In  any  analysis  of  LMS  there  is  the  problem  of 
nonstationarity  of  the  data  over  the  .entire  image  frame  that 
needs  to  be  considered.  This  problem  has  been  dealt  with  from 
a  variety  of  view  points.  In  addition  image  processing  has 
inherent  problems  caused  by  the  overlapping  of  data  blocks 
which  was  assumed.  In  the  original  derivation  of  the  one 
dimensional  LMS  algorithm,  independence  of  successive  input 
data  blocks  is  freguently  assumed.  This  is  not  necessarily 
true  for  image  data.  This  problem  has  been  addressed  in  the 
literature  under  the  assumption  of  specific  dependencies. 
These  range  from  independent  input  data  to  input  data 
possessing  considerable  correlation.  [Refs.  13-17]  An 
expansion  on  the  analysis  is  beyond  the  scope  of  this  thesis. 
Experimental  results  indicate  that  Widrow's  method  of  dealing 
with  nonstationarity  [Ref .  11] ,  as  adapted  by  Bitmead  [Ref . 
18]  and  Hsia  [Ref.  19]  and  termed  "The  Normalized  LMS  (NLMS),11 
has  some  general  applicability  in  the  processing  of  two- 
dimensional  images.  As  originally  presented  the  algorithm  was 
based  on  the  assumption  of  gaussian  independent  identically 
distributed  (iid)  data.  In  order  to  present  an  overview  of 
this  method,  as  adapted  to  two-dimensional  processes,  define 
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x(i)  = 


X(n~l,m) 

x(n-2 ,m) 


x(n-P+l,m) 


x(n,m-Q+l) 
x(n-l,m-Q+l) 


x(n-P+l,m-Q+l) 

(3-8) 
This  vector  consists  of  the  inputs  under  the  PxQ  prediction 
mask  at  update  i.   Further  define 


h(l,0) 
h(2,0) 


h(i)  = 


h(P-l,0) 


h(0 
h(l 


0-1) 
Q-l) 


h(P-l,Q-l) 

(3-9) 
as  the  prediction  coefficients  at  update  i  concatenated  into 
a  column  vector.  Using  the  vector  form  of  (3-7)  it  was 
proposed  by  Hadhound  and  Thomas  [Ref.  12]  that  the  actual 
coefficient  adaption  process  in  a  image  could  be  represented 
as 
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h°(i+l)  =  h°(i)  +  i(i) 

(3-10) 
The  term  2(i)  ^s  a  column  vector  of  random  increments  to  the 
ideal  time  varying  filter  which  would  track  the  input  data. 
Subtracting  (3-10)  from  (3-7)  in  vector  form  one  can  define 

V(i)  =  h(i)  -  h°(i) 

(3-11) 
which  represents  the  coefficient  deviation  from  the  ideal 
coefficients.  The  vector  V(i)  is  called  the  coefficient  error 
vector.  The  coefficient  error  update  can  be  written  as 
follows 

V(i+1)  =  (I-2/ix(i)x(i)T)V(i)  +  2/ix(i)e°(i)-l(i) 

(3-12) 
where  e°(i)  is  the  error  resulting  from  prediction  with  the 
optimum  time  varying  coefficients  h°(i).  The  coefficient 
vector  h°(i)  represents  the  optimal  coefficients  at  update  i 
as  determined  by  the  solution  of  the  normal  equations.  The 
ultimate  desire  is  for  the  coefficient  error  term  to  approach 
zero  as  i  becomes  large.  If  the  input  data  is  independent  and 
stationary,  this  can  be  accomplished  [Ref.  11]  with  a  step 
size  0<u<  A    ,  where  A    is  the  largest  eigenvalue  of  R. 

Am  ax  max 

However  in  the  current  form  (3-12)  it  can  be  seen  that  this 
would  not  produce  zero  coefficient  error  as  time  progressed. 
The  result  is  a  contribution  to  the  coefficient  deviation  due 
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to  the  approximation  of  the  gradient  (the  first  term  of  (3- 
12))  and  a  deviation  due  to  the  lag  in  tracking  the  ideal 
coefficient  changes  (the  last  terms  in  (3-12))  [Refs.  15  and 
16]  . 

In  the  "Normalized  LMS  (NLMS) "  algorithm  the  step  size  is 
defined  as 

H  =   _* — 
4Tr(R) 

(3-13) 

where  Tr(R)  may  be  approximated  by  the  ||x(i)  ||  .   The  form 

given  in  (3-12)  after  substituting  (3-13)  has  both  heuristic 

and  practical  appeal.   Both  terms  containing  the  input  data 

under  the  mask  have  been  normalized.  This  decreases  the 

explicit  dependence  of  the  coefficient  error  on  the  change  of 

data  input  power.   The  term  a  has  been  addressed  in  [Ref.  14] 

for  the  iid  case  and  related  to  the  ratio  of  the  variance  of 

the  step  change  in  (3-10)  and  the  MMSE.   This  is  of  interest 

from  a  theoretical  standpoint;  however,  in  practice  these 

factors  are  rarely  available.    In  practice  it  has  been 

observed  that  choosing  a  to  satisfy  0<a<2  for  the  case  of 

independent  data  works  well  with  the  most  rapid  convergence 

of  the  MSE  at  a=l.   In  Ref.  14  an  analysis  was  performed  on 

correlated  data  for  a  two  input  predictor.    The  results 

extended  to  the  "Normalized  LMS"  algorithm  would  indicate  that 

the  greater  the  correlation  of  the  input  data  the  smaller 
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should  be  the  a  term  (typically  a<l) .  After  substitution  of 
the  step  size  (3-13)  the  coefficient  update  in  vector  form 
becomes 

ae  ( n ,  m ;  i )  x  ( i ) 


h(i+l)  =  h(i)  + 

2Tr(R) 

(3-14) 

The  coefficient  error  results  in  excessive  MSE.   A  measure  of 

this  is  the  misadjustment  ratio  M  defined  as 

6(00)  -  £  . 

'      m  i  n 

M  =   


€     . 

m  l  n 

(3-15) 
where  e  (°o)  is  the  MSE  at  steady  state  and  e  .  is  the 
theoretical  MMSE  [Ref.  14]. 

Figure  3-2  shows  a  plot  of  the  misadjustment  ratio  versus 
various  values  of  the  a  term  in  the  step  size.  The  results 
presented  are  the  average  of  simulations  on  synthetic  images. 
The  image  process  is  modeled  by  a  2x2  quarter  plane  separable 
autoregressive  filter. 
x(n,m)  =  0.5  x(n-l,m)  -  0.4  x(n,m-l)  +  0.8  x(n-l,m-l)  +  w(n,m) 

(3-16) 

2 

The  results  are  for  a  range  of  a     from  0.6  to  4.0.   It  can  be 

seen  that  this  form  of  the  step  size  produces  a  unique  minimum 
M  for  a  given  process  model.  Considering  (3-12)  with  Figure 
3-2,  it  can  also  be  seen  that  as  a  becomes  small,  the  excess 
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MSE   increases  due  to  coefficient   lag.     Further,   as  a 
increases,  so  does  the  excess  MSE  due  to  gradient  error. 


Figure  3-2.   Misadjustment  Versus  a 


An  example  of  the  predictor  coefficient  deviation  is  shown 

in  Figure  3-3.   The  image  frame  processed  with  the  NLMS  PEF 

o 

was  synthesized  with  (3-16)  and  a  =0.64. 

The  initial  coefficients  were  set  by  a  least  squares 
technique  discussed  in  the  implementation  section.  It  can  be 
seen  for  the  processed  row  that  the  coefficients,  while  not 
the  same  as  those  for  the  original  process,  are  stable  and 
tracking.  The  row  processed  for  Figure  3-3  is  shown  in  Figure 
3-4. 
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Figure  3-3.   Predictor  Coefficient  Plot 
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Figure  3L4 

Process 


Row   of   the   Background 


The  desired  result  is  the  decorrelation  of  the  image 
process  (i.e.,  a  nearly  white  noise  error  residual);  this  can 
be  seen  in  Figure  3-5. 

Figure  3-6  contains  a  plot  of  the  normalized  histogram  of 
the  error  residual  output  of  the  NLMS  Prediction  Error  Filter. 
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Figure  3-5. 
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Error  Residue  Correlation 
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It  can  be  seen  that  the  error  residual  distribution  has 
a  gaussian  appearance  as  expected  with  a  gaussian  driving 
process. 

D.   LMS  PEF  IMPLEMENTATION 

One  of  the  deficiencies  in  the  LMS  algorithm  is  the  speed 

of  convergence.  The  detection  algorithm  discussed  in  a  later 
chapter  depends  on  the  ability  of  the  PEF  to  decorrelate  the 
background  process. 

Another  problem  that  is  inherent  in  the  processing  of  the 
image  is  that  of  directionality.  In  order  to  process  the 
image  a  scan  pattern  and  a  method  of  handling  boundary 
conditions  must  be  chosen.  Once  these  have  been  established, 
discontinuities  and  directional  dependencies  have  been 
introduced  into  the  process.  Although  directional  dependence 
cannot  be  overcome,  the  effect  of  discontinuities  and  settling 
time  can  be  decreased. 

The  processing  of  images  for  this  thesis  was  done  using 
a  PEF  with  a  3x3  rectangular  support  and  a  row  raster  scan 
pattern.  Initial  conditions  were  divided  into  two  parts.  The 
filter  coefficients  at  the  beginning  of  each  frame  were 
initialized  with  coefficients  derived  from  a  least  squares 
solution  for  a  small  area  (9x9)  pixels  (see  Figure  3-7)  at  the 
beginning  of  the  frame.  The  initial  coefficients  could  have 
been  derived  from  other  methods.  However,   the  method  used 


30 


c 
p 

o 

N 
a 

c 
0 

"~   U 

least   squares 
[OiiLQliz_Citjori_b_ox 

/ 

\ 

r 

0 

• 

• 

-*  n    .    .    .    . 

v 

direct 

on 

.    .    .    v  .    . 

/     5CQtl 

Figure  3-7.   Process  Scan  Diagram  With 
Initialization  Window 


should  have  a  misadjustment  less  than  or  equal  to  that 
anticipated  for  the  LMS  algorithm.  Initial  coefficients  for 
successive  rows  were  steepest  descent  algorithm  updates  of  the 
previous  row  initial  coefficients. 

The  remainder  of  the  frame  processing  follows  from  the 
solution  of  (3-2)  and  the  coefficient  update  computed  using 
(3-7).  The  online  computation  of  an  estimate  of  the  mis- 
adjustment  ratio  can  be  used  as  a  measure  of  the  algorithm 
efficiency.   Anomalous  jumps  in  the  misadjustment  can  be  used 
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to  signal  that  the  coefficients  should  be  reinitialized  using 
the  least  squares  method. 
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IV.     RECURSIVE  LEAST  SQUARES  FOR  TVO-DIMENSIONAL 

PREDICTION 

A.   THEORY 

The  least  squares  algorithm  is  a  method  used  in  the 

estimation  of  model  parameters.   [Refs.  6,  22,  and  23]   The 

objective  is  the  estimation  of  the  model  parameters  that  will 

fit  the  model  to  the  observed  data.    The  criterion  for 

determining  the  goodness  of  fit  is  the  minimum  sum  of  the 

squared  error  between  the  model  output  and  all  observations. 

An  extension  of  this  method  is  the  Recursive  Least  Squares 

(RLS)   algorithm.    This  method  is  identical  to  the  least 

squares  algorithm  but  provides  for  updating  the  estimate  of 

the  model  parameters  as  more  data  becomes  available.    In 

Chapter  II  the  concept  of  fitting  a  reduced  order  model 

(prediction  model)  to  an  image  was  addressed.  This  model  will 

be  used  in  the  RLS  algorithm  and  has  the  form 

x(n,m)  =  E  S  h(£,k)  x(n-£,m-k)  +  e(n,m) 
£  k 

(£,k)  e   a 
(£,k)  f    (0,0) 

(4-1) 
The  term  x(n,m)  is  the  observation  and  the  terms  x(n-£,m-k) 
are  the  past  data.  The  terms  h(£,k)  are  the  unknown 
parameters  that  must  be  determined  in  order  to  fit  the  model 
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to  the  observation.   The  term  e(n,m)  is  the  error  between  the 

model  and  the  observed  image  data  at  position  (n,m)  .   The 

least  squares  criterion  is  as  follows 

N-l  M-l 
J(h,x;N,M)  =   S   S   eZ(h,x;n,m) 
n=0  m=0 

(4-2) 
where  J(h,x;N,M)  is  the  cost  function  to  be  minimized  by  an 
appropriate  choice  of  h(£,k;N,M).  The  error  term  in  (4-2)  is 
defined  as 

e(h,x;n,m)  =  x(n,m)  -  x(n,m) 

(4-3) 

The  term  x(n,m)  is  the  prediction  of  the  point  x(n,m)  defined 

as 

P-l  Q-l 
x(n,m)  =   2    S   h(£,k;N,M)  x(n-£,m-k) 
£=0  k=0 
(ffk)*(ofo) 

(4-4) 
where  the  support  of  the  predictor  model  is  PxQ.  The  second 
indices  (N,M)  in  the  term  h(£,k;N,M)  are  included  to  indicate 
that  h(£,k;N,M)  is  held  constant  for  a  block  of  data  0<n<N- 
1  and  0<m<M-l.  Thus  it  is  the  sum  of  the  errors  with  this 
choice  for  the  coefficients  that  is  minimized. 

The  criterion  as  defined  above  in  conjunction  with  the 
prediction  model  can  be  recast  in  the  form  of  a  PEF.  Figure 
4-1  shows  a  3x3  rectangular  PEF  as  defined  by  (4-5) . 
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Figure  4-1.   Prediction  Error  Filter 


P-l  Q-l 
e(n,m)  =  x(n,m)      E   S   h(£,k;N,M)  x(n-£,m-k) 

£=0  k=0 
(trk)*(0f0) 

(4-5) 
The  output  of  the  PEF  is  the  error  term  to  be  minimized  as 
defined  by  the  cost  function.  In  order  to  simplify  the 
notation  the  h  and  x  variables  have  been  dropped  from  the 
error  term;  however,  the  dependency  is  understood.  As 
proposed  in  Chapter  II,  if  the  error  term  is  minimized  and  the 
original  process  was  a  gaussian  WNDR,  then  the  error  will  be 
white  and  gaussian  distributed. 
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At  this  point  it  is  advantageous  to  define  the  previous 
expressions  in  vector  form.  Although  the  ultimate  objective 
and  results  remain  the  same,  directionality  is  inherent  in  the 
form  to  follow.  The  development  to  follow  will  be  for  a 
raster  scan  by  row,  but  it  may  be  generalized  for  any  scan 
pattern  for  which  recursive  computability  is  assured.  The 
mapping  function  for  a  row  scan  (see  Figure  4-1)  is 

i  =  n  +  Nm   0<n<N-l,  0<n<N-l 

(4-6) 
where  N  represents  the  row  length  and  M  represents  the  column 
length  of  the  data  field  defined  on  a  rectangular  lattice. 
By  using  this  mapping  function  to  replace  the  indices  (n,m) 
and  concatenating  terms  by  rows  into  a  column  vector  the 
following  terms  can  be  redefined. 


x(i)  = 


x(n-l,m) 
x(n-2 ,m) 


x(n-P+l,m) 


x(n,m-Q+l) 
x(n-l,m-Q+l) 


x(n-P+l,m-Q+l) 


and 


(4-7) 
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h(j)  = 


h(l,0) 
h(2,0) 


h(P-l,0) 


h(0,Q-l) 
h(l,Q-l) 


h(P-l,Q-l) 

(4-8) 
The  vector  x(i)  is  thus  comprised  of  the  points  under  the 
prediction  mask  at  update  i=n+Nm.  The  vector  h(j)  is 
comprised  of  the  prediction  filter  coefficients  based  on  all 
data  for  0<i<j-l.  It  follows  from  this  mapping  that 
x(i)  =  x(n,m)  and  x(i)  =  x(n,m) 


(4-9) 


and 


x(i)  =  hT(j)  x(i) 


0<i<j-l 


(4-10) 


The  cost  function  is  then  expressed  as 

j-1  j-(i+l)    2   j-1  j-(i+l) 


J  ( h ,  x ;  j  )  =  2  X 
i=0 


e(i)c  =  2  X 
i=0 


(x(i)-h'(j)x(i)) 


(4-11) 
Where  the  term  xj"(l+1)  is  added  to  give  the  cost  function  a 
fading  memory  [Ref.  23]  and  thus  a  greater  ability  to  adapt 
to  nonstationarities  in  the  image.   In  a  stationary  image  the 
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exponential  weighting  factor  would  be  set  to  unity.  From  this 

point  on,  the  development  of  the  RLS  algorithm  will  parallel 

Haykin's  one  dimensional  approach  in  Refs.  22  and  23. 

The  cost  function  as  presented  is  quadratic  in  the  terms 

of  h(j)  and  therefore  possesses  a  unique  minimum.   To  solve 

for  the  optimal  coefficients  h(j)   in  terms  of  the  cost 

function,  (4-11)  is  differentiated  with  respect  to  the  filter 

coefficients. 

j-1  j-(i+l)       j-1  j-(i+l) 
3J(h,x;:j)  =  _2   2  X  x(i)  -     2  X  h(j)x(i) 

3h(j)        i=0  i=0 

(4-12) 

Now  as  in  Ref.  23  the  following  definitions  will  be  made. 

The  deterministic  correlation  matrix  is  defined  as 

j-1  j-(i+l) 
*(j)  =  2  X        x(i)x(i) 
i=0 

(4-13) 

In  the  same  vein  the  deterministic  cross  correlation  vector 

(i.e. ,  the  cross  correlation  of  the  point  to  be  predicted  with 

the  data  vector)  is  defined  as 

j-1  j-(i+l) 
£(j)  =  2  A        x(i)x(i) 
i=0 

(4-14) 

Using  the  last  two  definitions,  setting  the  derivative  of  the 

cost  function   (4-12)  equal   to  zero  and   rearranging  yields 
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*(j)  h(j)  =  0(j) 

(4-15) 
This  form  has  been  termed  in  Refs.  22  and  23  as  the 
deterministic  normal  equations.  The  solution  of  (4-15) 
produces  the  optimal  coefficients  h(j)  in  the  least  squares 
sense. 

h(j)  =  *_1(j)  *(j) 

(4-16) 

The  objective  now  is  to  solve  for  h(j)  recursively.   This 

can  be  accomplished  by  forming  the  right  hand  side  of  (4-16) 

recursively.   [Ref.  22]   From  (4-13)  the  correlation  matrix 

can  be  expressed  in  the  recursive  form 

*(j)  =  A*(j-1)  +  x(j)x(j)T 

(4-17) 
Further,  from  (4-14) ,  the  cross  correlation  vector  may  be 
formed  recursively  as 

£(j)  =  A  0(j-l)  +  x(j)x(j) 

(4-18) 
In  forming  the  inverse  of  the  correlation  matrix  the  matrix 
inversion  lemma  [Ref.  22]  is  used.  On  application  of  the 
lemma,  the  inverse  of  the  correlation  matrix  becomes 


*-!(j)  =  VV^j-l)  - 


V2#  ^j-l)  x(j)x(j)T*  !(j-l) 
1  +  X_1xT(j)  f_1(j-l)x(j) 


(4-19) 
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From  (4-19)   two  terms  are  defined  which  will  facilitate 
computations.   These  are 

P(j)  =  *_1(j) 

(4-20) 

and 

A-1  P(j-l)  x(j) 


E(D)  = Tt 

1  +  X   V(j)  P(j-l)x(j) 

(4-21) 
where  K(j)  is  called  the  "gain  vector."  The  form  of  K(j)  can 
be  shown  [Ref.  22]  to  reduce  to 

K(j)  =  P(j)  x(j) 

(4-22) 
It  also  follows  from  (4-19,  20,  21  and  22)  that 

P(j)  =  A"1  P(j-l)  "  *-1K(J)  xT(j)  P(j-l) 

(4-23) 
The  update  of  the  filter  coefficient  vector  h(j)  can  now  be 
formed  with  the  use  of  (4-16,  17,  21,  and  22). 

h(j)  =  h(j-l)  +  K(j)  [x(j)  -  hT(j-l)  x(j)] 

(4-24) 
Using  this  equation  the  prediction  error  is  defined  as 

e(j)  =  x(j)  -  hT(j-l)x(j) 

(4-25) 
This  has  been  termed  the  a  priori  prediction  error.  Further 
define  the  a  posteriori  prediction  error  as 
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a(j)  =  x(j)  -  hT(j)x(j) 

(4-26) 
It  is  the  a  prior  prediction  error  that  is  used  in  the  cost 
function  in  (4-2)  .   Now  using  (4-25)  in  (4-24)  the  final  form 
of  the  coefficient  update  equation  is 

h(j)  =  h(j-l)  +  K(j)e(j) 

(4-27) 
This  equation  completes  the  necessary  equations  for  the 
construction  of  the  adaptively  updated  prediction  error 
filter. 

A  term  that  will  be  useful  in  the  detection  algorithm 
(covered  in  the  next  chapter)  is  the  minimum  value  of  the  cost 
function.  The  cost  function  attains  the  minimum  value  if  at 
each  step  in  the  recursion  the  filter  coefficients  satisfy  the 
deterministic  normal  equations  (4-15) .  In  this  case  the  cost 
function  may  be  written  as 

J  ■  (J)  =  V"(i+1)  J  .  (j-1)  +  e(j)a(j) 

(4-28) 
The  weighted  average  of  Jmin(j)  over  the  range  of  the 
iterations  modified  by  the  exponential  weight  factor  provides 
an  estimate  of  the  error  variance  at  step  j .  The  effective 
memory  length  for  averaging  purposes  is 

MEM  =  _*- 
1-X 

(4-29) 
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B.  ALGORITHM  ANALYSIS 

The  properties  of  the  least  squares  method  of  parameter 

estimation  have  been  addressed  extensively  in  Ref .  23  and  Ref . 
24.  In  Chapter  II  it  was  assumed  that  the  background  process 
could  be  modelled  as  a  WNDR.  It  was  further  assumed  that,  if 
our  reduced  order  model  was  sufficiently  large  and  a  minimum 
variance  fitting  technique  was  employed,  the  resultant  error 
would  be  white.  Under  these  assumptions  the  least  squares 
estimate  of  h(£,k;M,N)  has  some  important  properties,  namely 
that  it  is  consistent  and  the  best  linear  unbiased  estimator. 
This  has  been  shown  for  the  one-dimensional  least  squares 
problem  in  [Refs.  23,  24,  and  25]  and  extends  readily  to  the 
two-dimensional  problem.  Further  it  has  been  shown  [Ref.  2  2 
and  23]  for  the  one  dimensional  case  that  the  mean  square 
error  will  converge  in  approximately  twice  the  filter  length 
to  the  minimum  mean  square  error.  This  was  shown 
experimentally  to  be  a  good  estimate  for  the  two  dimensional 
case  as  well. 

In  addition,  if  the  driving  source  for  the  model  is 
gaussian,  then  the  least  squares  estimate  results  in  the 
maximum  likelihood  estimate.  In  order  to  demonstrate  this 
some  terms  will  need  to  be  defined.  Assume  as  before  that 
the  error  residual  is  white  and  gaussian.  The  joint 
probability  distribution  for  the  error  residual  in  (4-11) , 
without  the  exponential  weighting  factor  is 
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p(e(0) ,e(l) ,...e(j-l))  = 


T~i i/2   exp 


j_1 

-2  e(i)2 

i=o 
~~2 


20 


(4-30) 
Inserting  (4-25)  for  the  term  e(i)  the  likelihood  function  can 
be  written  as 


Lx(h(j),a/)  = 


^7"  i/2  exp 


j_1         T        2 
2  (x(i)  -  h'(j)x(i))^ 


-  i=o 


20 


(4-31) 
It  can  be  seen  that  maximizing  the  likelihood  of  x(i),  the 
observations,  with  respect  to  the  parameter  h  to  be  estimated 
is  the  same  as  the  least  squares  criteria. 

The  algorithm  has  its  faults  as  well  as  its  good  points. 
Without  the  exponential  weighting  factor  the  memory  length  is 
unlimited.  This  presents  a  problem  in  implementation  on  fixed 
word  length  devices.  Normalization  of  the  P  matrix  has  been 
attempted  to  overcome  this  problem;  however,  this  introduces 
a  bias  in  the  mean  square  error.  The  long  memory  also  makes 
the  algorithm  less  responsive  to  changes  in  the  image  as  time 
progresses.  This  was  the  reason  for  the  inclusion  of  the 
exponential  weighting  factor  in  the  recursions.  There  is  a 
drawback  to  the  use  of  the  weighting  factor  however.  The 
smaller  the  factor  (i.e. ,  shorter  the  memory) ,  the  greater  the 
error  variance.  Experimental  results  with  synthetic  images 
have  shown  this  to  be  true.   At  a  memory  length  of  less  than 
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approximately  ten  pixels  the  results  were  indistinguishable 
from  those  of  LMS. 

C.   RLS  PEF  IMPLEMENTATION 

The  following  is  the  complete  RLS  algorithm.   Note  that 

the  initialization  given  here  is  the  traditional  one.  The 
actual  initialization  techniques  employed  in  this  thesis  are 
covered  following  the  algorithm  listing.  The  initial 
conditions  on  P.  h,  and  J  •  are 

'  — '        mm 

P(0)  =  6_1I 
h(0)  =  0 
J(0)  =  0 

(4-32) 
where  <S  is  a  small  number  (approx.  <<1)  .  The  exponential 
weight  factor  is  set  at  the  desired  value  in  the  range  0<X<1. 
The  update  recursion  for  i>0  is  then 

A-1  P(i-l)x(i) 


-(1)  =  — n 

1  +  A  V  (i)P(i-l)x(i) 

(4-33) 
where  K(i)  is  the  gain  vector.  The  inverse  correlation  matrix 
P  is  updated  by 

P(i)  =  A_1P(i-l)  -  A_1K(i)xT(i)P(i-l) 

(4-34) 
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The  a  priori  prediction  error  is 


e(i)  =  x(i)  -  hT(i-l)x(i) 


(4-35) 


(4-36) 


(4-37) 


The  filter  coefficients  are  updated  as  follows 

h(j)  =  h(j-l)  +  K(j)e(j) 

The  posteriori  prediction  error  update  is 

a(i)  =  x(i)  -  hT(j)x(i) 

The   final    item  to  be  computed   is  JBl-n 

J    .    (i)    =  \J   .    (i-1)    +   e(i)a(i) 

min  min  * 

(4-38) 

It  is  from  this  last  term  that  an  estimate  is  formed  of  the 

error  variance.   This  estimate  is 

o  2   =   J  .  (i)/(n-8) 

e       m  i  n  x  '  '    x     ' 

(4-39) 

where  n  is  the  number  of  points  processed  and  the  n-8  term  is 
to  account  for  the  loss  of  8  degrees  of  freedom  due  to 
estimating  the  filter  coefficients  from  the  data  set  [Ref. 
24]. 

The  speed  of  convergence  has  been  previously  mentioned. 
While  a  convergence  rate  of  twice  the  filter  length  may 
provide  an  acceptable  rate  of  convergence  for  continuous  data 
sets  it  presents  a  problem  in  image  processing.  The  problem 
arises  in  three  areas: 
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1)  at  the  initial  start  up  on  each  frame, 

2)  at  the  beginning  of  each  row  or  column  in  a  row  or 
column  raster  scan,  and 

3)  where  a  large  change  in  background  statistics  occurs. 
In  order  to  reduce  the  effects  of  the  problem  areas  the 

least  square  algorithm  was  implemented  in  the  form  of  a 
secondary  processor.  Logic  flags  were  established  in  the  RLS 
processing  to  signal  the  occurrence  of  one  of  the  problems. 
Once  one  of  the  flags  is  raised  the  RLS  algorithm  is  re- 
initialized with  the  least  squares  algorithm.  The  least 
squares  algorithm  is  applied  to  a  9x9  rectangular  block  of 
data  centered  on  the  flagged  point  in  Figure  4-2.  The 
algorithm  produces  an  estimate  of  the  filter  coefficients,  the 
P  matrix,  and  Jmin  which  are  passed  to  the  main  RLS  processing 
routine. 

It  can  be  seen  that  this  procedure  can  introduce  a 
considerable  burden  in  the  processing.  A  suitable  alternative 
is  to  eliminate  the  row  initialization  flag  while  still 
retaining  the  frame  initialization  flag.  The  initial  estimate 
of  P,  h(j),  and  Jmin  for  each  row  or  column  could  be  obtained 
from  the  initial  value  of  the  previous  row  or  column.  The 
update  is  repeated  for  successive  rows  or  columns  in  the  same 
manner. 
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Figure  4-2.   RLS  PEF  Initialization 

An  example  of  the  experimental  results  follows  using  the 
two-dimensional  RLS  algorithm.    The  image  processed  was 
synthetically  produced  using  a  2x2  autoregression  mask.   The 
actual  process  model  was  separable  and  defined  as 
x(n,m)  =  0.5x(n-l,m)  -  0.4x(n,m-l)  +  0. 8x(n-l,m-l)  +  v(n,m) 

(4-40) 

2  2 

with  the  variance  of  the  source  a  =2 . 0  and  output  o   =14.1. 

V  X 

The  overall  mean  of  the  image  was  removed  prior  to  processing. 
The  correlation  function  of  the  error  residual  resulting  from 
processing  the  entire  image  is  shown  in  Figure  4-3. 
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Figure  4-3.   Error  Residual  Correlation 
Function 


Figure  4-4  is  a  histogram  plot  of  the  error  residual  with 
a  gaussian  distribution  superimposed  for  comparison. 

Figure  4-5  is  a  plot  of  the  filter  coefficients  resulting 
from  a  RLS  row  raster  scan  of  the  same  row  as  in  Figure  4-3. 

This  concludes  the  development  and  implementation  of  the 
two  dimensional  RLS  algorithm.  The  RLS  algorithm  was  found 
to  be  superior  to  the  LMS  algorithm  in  all  areas  except  in  the 
area  of  computations  required.  The  performance  will  prove  to 
be  an  overriding  consideration  in  the  target  detection 
problem. 
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Figure  4-4.   Histogram  of  Error  Residue 
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Figure  4-5.   Predictor  Coefficient  Plot 
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V.   TARGET  DETECTION 

A.   GENERAL 

The  previous  three  chapters  have  dealt  with  the  theory  and 
algorithms  for  modelling  the  background  process.  It  was  shown 
that  the  result  of  inverse  filtering  with  an  adaptive  PEF  is 
a  nearly  white  error  residual  process.  It  is  this  error 
residual  that  is  of  particular  interest  in  the  detection 
algorithm.  The  residual  represents  the  unpredictable  portion 
of  the  image  process.  The  question  that  arises  is  the  effect 
of  the  resultant  decorrelation  on  point  targets  present  in  the 
image. 

When  a  target  is  present  in  an  image,  it  replaces  the 
background  pixel  at  that  location.  Therefore  the  target  is 
opaque  as  far  as  the  image  process  is  concerned.  Further  the 
target,  if  present,  is  deterministic  and  its  intensity  remains 
constant  over  time.  The  probability  density  of  the  image 
versus  the  target  may  be  viewed  as  in  Figure  5-1.  In  this 
figure  the  intensity  of  the  target  is  a  fixed  value  T  and  its 
PDF  is  represented  as  an  impulse;  the  mean  of  the  PDF  for  the 
background  is  specified  by  m.  It  can  be  seen  that  only  the 
magnitude  of  the  difference  between  the  image  mean  and  the 
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Figure   5-1.     Background   PDF   Versus 
Imbedded  Target 


target  intensity  T  is  notable  when  compared  to  the  background 
process.  Therefore  prior  to  processing  with  the  adaptive  PEF 
the  background  process  mean  is  removed.  This  results  in  a  new 
target  T  with  an  intensity  of  Am. 

During  the  modelling  phase  the  coefficients  of  the  reduced 
order  model  are  adapted  to  fit  the  image  at  each  pixel 
location.  The  error  in  the  fit  is  relatively  small  in  the 
locations  occupied  by  actual  background  pixels.  This  also 
holds  true  for  any  pixel  value  that  closely  resembles  a 
background  pixel;  a  target  with  a  small  Am  is  more  difficult 
to  distinguish  from  the  background  than  one  with  a  larger  Am. 
The  value  of  the  error  process  at  the  target  location  is 
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P-l  Q-l 

e*(n,m)  =  T*(n,m)  -22     h(£,k)x(n-€,m-k) 

i=0   k=0 
(«,k)*(0,0) 

(5-1) 

The  term  e  (n,m)  represents  the  error  due  to  the  presence  of 

the  target.   This  term  may  be  rewritten  as  a  combination  of 

two  error  terms.   The  first  term  represents  the  error  that 

would  be  result  from  predicting  the  background  in  the  absence 

of  a  target.   That  is, 

P-l  Q-l 
e(n,m)  =  x(n,m)  -  2   2  h(£,k)  x  (n-£,m-k) 

£=0  k=0 

(5-2) 
The  other  term  is  the  difference  between  the  target  pixel  and 
the  missing  background  pixel. 

e-r(n,m)  =  e  (n,m)  -  e(n,m)  =  T  (n,m)  -  x(n,m) 

(5-3) 
The  resultant  error  process  e  (n,m) ,  separated  in  this  manner, 
can  be  compared  to  the  background  error  process.  The  process 
e(n,m) ,  assuming  complete  decorrelation  of  the  background,  can 
be  completely  characterized  by  its  variance.  On  the  other 
hand  the  error  process  e-r(n,m)  possesses  a  mean  which  will  be 
non-zero  and  is  given  by 

E{e-r(n,m)  }=  Am  =  T*(n,m) 

(5-4) 


53 


The  variance  of  error  process  e  (n,m)  is  defined  as 
Var{e*(n,m)}  =  E{e*(n,m)2}  -  E{e*(n,m)}2 

(5-5) 
After  using  (5-1)  through  (5-4)  in  (5-5)  and  simplifying  the 

variance  of  the  target  error  process  becomes 

*  2  2 

Var{e  (n,m) }  =  a^    (n,m)  -  ag  (n,m) 

(5-6) 

Intuitively  this  means  that  the  ability  to  recognize  the 

target  is  linked  to  the  contrast  between  the  background  random 

process  and  the  deterministic  target.    These  two  error 

processes  present  in  a  prediction  error  image  may  be  viewed 

as  in  Figure  5-2.   It  can  be  seen  that  the  magnitude  of  the 

mean  difference  between  the  target  pixel  and  the  anticipated 

background  pixel  determines  the  separation  between  the  two 

density  functions.   It  must  be  kept  in  mind  that  the  density 

function  is  for  the  local  area  since  the  filter  is  spatially 

adaptive.   Thus  it  is  the  difference  between  the  local  area 

mean  and  the  target  intensity  that  constitutes   (Am)  the 

difference  between  the  two  density  functions.   The  factor  is 

evident  in  an  image  with  multiple  background  processes.   The 

solution  presented  was  for  a  homogenous  (single)  background 

process.   In  an  image  with  multiple  background  processes  the 


Note  that  as  previously  discussed  x(n,m)  and  therefore 
e(n,m)  has  a  mean  of  zero. 
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Figure  5-2 .    Background  Error  Process 
Versus  Target  Error  Process 


removal  of  the  image  mean  does  not  necessarily  produce  a  local 
process  PDF  with  a  mean  of  zero.  A  separation  between  the 
local  background  process  and  the  target  process  still  results 
but  it  is  now  with  reference  to  the  local  process  mean. 

The  form  of  the  inverse  filter  as  mentioned  previously  is 
a  prediction  error  filter  (PEF) .  This  filter  has  a  finite 
impulse  response.  It  follows  that  the  effect  of  a  target  on 
the  error  process  is  not  limited  to  the  pixel  where  the  target 
is  located.  The  PEF  used  in  this  thesis  is  rectangular  with 
a  finite  impulse  response  of  3x3  pixels  in  dimension.  If  the 
background  process  is  stationary  within  this  region,  the  error 
process  can  be  defined  as  follows.  The  error  at  each  location 
within  the  impulse  response  support  region  is 
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2         2 

e*(n+i,m+j)    =     22     a (£,k;n+i,m+j )    x(n+i-£ ,m+j-k) 
£=0   k=0 

<«,k)*(i,j) 

+   a(i,m;n+i,m+j)T   (n,m) 

(5-7) 
where  a(£ , k;n+i,m+j )    are  the   filter  coefficients 

a(«,k;n+i,m+j)    =1  1  (^k)    =    (0,0) 

-h(£,k)  0<£<2,    0<k<2 

(5-8) 
The  additional  indices  (n+i,m+j)  are  to  indicate  that  these 
coefficients  are  spatially  variant.  Using  (5-3) ,  (5-7) ,  and 
(5-8)  and  concatenating  the  respective  components  into  column 
vectors,  the  error  process  in  the  support  region  may  be 
written  as 

e  =  e  +  e-r(n,m)a 

(5-9) 
Now  assume  that  the  filter  is  at  steady  state  within  the 
region.  Since  the  background  has  mean  zero,  e  has  mean  zero. 
Then  the  mean  of  the  error  process  at  each  location  within  the 
region  in  vector  form  follows  from  (5-4)  and  (5-9) 

E{e  )  =  Ama 

(5-10) 
Since  the  error  process  is  uncorrelated,  its  covariance  matrix 
K  is  diagonal  and  follows  from  (5-6)  and  (5-8)  .  The  diagonal 
elements,   k(i,j),  of  the  covariance  matrix  are  defined  as 
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r 


*e(i,j)  =  < 


2     2  •  . 

x       e 

a(i,j)2<7x2  +  oeZ  i=J^O 


(5-11) 

The  purpose  for  the  preceding  derivation  of  the  error 
process  in  the  impulse  region  is  twofold.  First,  it  provides 
some  additional  information  which  could  be  used  for  methods 
of  detection  more  oriented  to  hypothesis  testing  or 
significance  testing  of  multiple  pixel  samples.  Second,  it 
shows  that  some  significant  error  residual  will  occur  in  the 
immediate  area  of  the  target  in  the  impulse  response  support 
region.  The  error  is  the  result  of  the  adaptive  process  that 
attempts  to  minimize  the  error  variance.  It  seems  reasonable 
that  the  greatest  power  in  the  error  residual  in  the  impulse 
response  support  region  will  be  concentrated  in  the  pixels 
which  are  the  nearest  neighbors  of  the  target.  This  concept 
can  be  used  in  grouping  suspected  targets  (results  of  the 
significance  test)  and  for  track  file  reduction. 

B.   SIGNIFICANCE  LEVEL 

The  concept  of  significance  testing  was  discussed  in 

Chapter  II.  Recall  that  a  significance  test  is  based  on  a 
null  hypothesis  Hq  whose  distribution  is  known.  It  is  against 
this  hypothesis  that  all  other  events  are  tested.  The  null 
hypothesis  HQ  to  be  used  in  this  section  is  that  of  the  error 
residual  process.   If  the  background  process  is  gaussian  then 
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the  error  process  is  also  gaussian.   In  the  preceding  section 

the  error  process  was  discussed.  The  error  process  due  to  the 

background  was  defined  in  (5-2)  and  the  assumption  was  made 

that  the  process  e(n,m)  was  completely  uncorrelated.   The 

probability  density  function  for  this  process  can  thus  be 

written  as 

2 


p(e(n,m))  =  —± exp 


-1  e(n,m)  ,. 
_  2  CTe(n,m)^_ 


(5-12) 
The  term  a  (n,m)  is  indexed  to  reflect  its  dependence  on  the 
local  statistics  of  the  image.  This  density  function  can  be 
used  to  define  the  significance  level. 

Since  the  prediction  error  residual  represents  the 
unpredictable  portion  of  the  background  it  contains  all  of  the 
pertinent  information,  without  redundancy,  necessary  for 
making  decisions.  For  this  reason  significance  testing  will 
be  performed  on  the  prediction  error  residual.  The  ordered 
set  of  all  such  residuals  will  be  referred  to  as  the 
prediction  error  image. 

It  is  desirable  to  set  the  significance  level  in  such  a 
way  that  little  or  no  adjustment  is  needed  over  the  entire 
image.  To  start  with,  a  decision  is  made  to  reject  Hq  (i.e., 
to  consider  that  a  target  may  be  present)  whenever  the 
magnitude  of  the  error  falls  within  a  critical  region.  The 
lower  boundary  of  this  region  is  defined  by  the  maximum 
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acceptable  two  tailed  probability  of  a  type  I  error,  P  (reject 
Hq  given  Hq  is  true)  (see  Figure  5-3).  The  computed 
probability  of  a  type  I  error,  based  on  the  null  hypothesis, 
will  be  called  a.  Associated  with  this  probability  is  a  level 
B  which  defines  the  critical  region  (see  Figure  5-3)  . 


Figure   5-3.     Diagram  of  Significance 
Testing  Terminology 


Using  these  last  two  ideas  the  significance  level  a  may 
be  written  as 


00 

(27r)1/2  0    (n, 
*  B 


m) 


exp 


1  e(n,m)2 

2  a   (n,m) 

e  \    / 


de (n,m) 


(5-13) 
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and  the  corresponding  critical  region  is  defined  by 

1*1  >  B 

(5-14) 

It  can  be  seen  that  (5-13)  and  (5-14)  depend  on  the  local 

statistics.    If  a  fixed  threshold  B  (see  Figure  5-3)  is 

chosen,  the  probability  a  would  fluctuate  with  the  changes  in 

2  ...  ... 

variance  a    (n,m) .   This  is  undesirable.   The  probability  a  is 

equal  to  total  area  under  the  PDF  within  the  critical  region. 

The  goal  is  to  find  a  threshold  which  will  produce  a  constant 

a   independent  of  the  local  statistics. 

A  solution  is  to  map  the  error  distribution  into  a 

normalized  distribution  with  the  mapping 

z  =  e(n,m) 
ffe(n,m) 

(5-15) 

This  mapping  results  in  a  unit  variance  normal  distribution 

in  the  variable  z  and  the  critical  region  is  mapped  into  the 

region 


z I  >  A    where  A  =  — ^ 


ae    (n,m) 


(5-16) 


The  probability  a   is  determined  by 
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— 1 

(27T) 


exp 


:  iz2i 


1/2         r[      2 


iz  dz 


(5-17) 
A  depiction  of  the  mapped  error  residual  PDF  is  shown  in 
Figure  5-4. 


Figure   5-4.      Thresholding   in   the 
Normalized  PDF  Space 


Now  observe  that  a  choice  of  a  threshold  value  of  A=constant 

produces  a  constant  probability  a.   The  test  (5-16)  can  thus 

be  expressed,  using  (5-15)  as 

2 


eljun)    >  a2  =  constant 
o   (n,m) 


(5-18) 
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From  this  it  can  be  seen  that  a  constant  alarm  rate  (a)  can 
be  achieved  with  a  constant  threshold.  The  implication  of 
this  is  that  the  statistic  on  the  left  hand  side  of  (5-18)  is 
independent  of  the  variations  in  the  background  or  error 
residual  image.  Ideally  then,  a  single  threshold  would  be 
suitable  for  the  entire  image.  This  constant  threshold  may 
be  found  experimentally  by  adjusting  the  threshold  in  the 
prediction  error  images  to  produce  on  the  average  a  specified 
percentage  of  false  alarms  (pixels  which  exceed  the 
threshold) . 

The  preceding  discussion  defined  the  statistic  to  be 
tested  and  showed  that  the  threshold  may  be  set  to  any 
constant  which  produces  the  desired  constant  false  alarm  rate. 

The  terms  e(n,m)  are  computed  using  either  the  LMS  or  RLS 

2 

method.   In  addition  the  term  a    (n,m)  needs  to  be  computed. 

In  the  RLS  algorithm  this  estimate  was  formed  by  (4-39)  . 

2 

During  the  LMS  process  the  estimate  of  a    (n,m)  is  computed 

from  the  difference  in  the  unbiased  spatial  averages  of  the 

2  2 

input  x(n,m)   and  the  prediction  x(n,m)  .   To  compute  this  an 

unbiased  estimate  of  the  correlation  matrix  (R  )  was  used  in 

conjunction  with  a  column  vector  of  the  PEF  coefficients  as 

discussed  earlier.    The  following  equation  defines  the 

estimate  of  the  error  variance  and  is  computed  at  each  pixel 

location. 
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2  T 

o   (n,m)  =  a  R  a 

e  —    x  — 

(5-19) 

This  completes  the  description  of  the  significance  test. 

A  block  diagram  of  the  process  through  the  testing  stage  is 

shown  in  Figure  5-5. 
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Figure  5-5.   Diagram  of  Detection  Process 

The  testing  can  be  carried  out  at  the  same  time  the  image  is 
being  whitened.  This  reduces  the  requirements  for  additional 
storage.  It  was  mentioned  earlier  that  the  significance 
testing  should  be  followed  by  a  method  for  managing  the 
resultant  target  list.  This  will  not  be  addressed  in  this 
thesis.  However  a  possible  tracking  algorithm  is  covered  in 
the  next  chapter. 

C.   SIMULATION  RESULTS 

The  simulations  conducted  using  the  processing  algorithms 
developed  up  to  this  point  were  on  images  with  homogenous 
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backgrounds.   The  image  mean  was  computed  and  removed  prior 

to  processing;  however,  the  local  mean  was  not  removed.   The 

background  of  the  image  was  generated  from  one  of  several 

models  of  the  form 

2    2 
x(n,m)  =  2   2  h(i  ,k)  x(n-£,m-k)  +  w(n,m) 
£=0  k=0 
(«,k)*(0,0) 

(5-20) 

where  h(£,k)  represents  a  fixed  set  of  filter  coefficients  and 

w(n,m)  is  a  gaussian  white  noise  source  with  a  variety  of 

variances.  After  preparing  the  image  background,  targets  were 

placed  at  random  locations  within  the  image.  The  significance 

testing  results  with  the  various  models  and  driving  source 

noise  were  consistent  throughout.   The  dimensionality  of  the 

PEFs  used  were  3x3  pixels;  this  was  more  than  adequate  for  the 

2x2  model.   It  should  be  realized  that  in  other  cases  the  size 

of  the  PEF  required  to  whiten  a  specific  image  may  be  larger 

or  smaller  than  3x3. 

The  specific  results  presented  in  this  chapter  are  for  the 

background  model  in  (5-20)  with  the  coefficients  h( 1,0) =0.5, 

h (0,1) =0.8,  h( 1,1) =-0.4  and  a2  =2.0.   The  imbedded  targets 

were  set  to  specific  dB  levels  relative  to  the  local  variance 

in  the  9x9  pixel  region  in  which  they  were  to  be  imbedded. 

This  has  been  designated  TBR  (target  to  background  ratio) . 

Once  the  TBR  was  set,  the  mean  of  the  9x9  region  was  added  to 
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the  target.  The  target  then  replaced  the  pixel  at  the  center 
of  the  9x9  region.  The  purpose  of  this  was  to  test  the 
capability  of  the  whitening  and  detection  algorithms  on 
targets  that  differed  by  very  small  amounts  from  their 
neighbors.  Some  simulations  were  conducted  for  images  with 
imbedded  targets  corrupted  by  additive  white  noise.  The  level 
of  noise  is  specified  as  the  BNR  (background  to  noise  ratio) . 
The  averaged  results  for  multiple  simulations  at  each 
parameter  setting  are  shown  in  the  Appendix. 

Plots  of  the  results  of  processing  identical  images  with 
the  LMS  algorithm  and  the  RLS  algorithm  are  shown  in  Figures 
5-6  and  5-7  respectively.  These  plots  are  for  the  conditions 
of  no  noise  and  a  TBR  of  5  dB  with  the  threshold  set  to  11.4 
and  ten  targets  imbedded  in  the  image.  Figures  5-8  and  5-9 
are  the  same  images  with  imbedded  targets  of  the  same  TBR  as 
in  Figures  5-6  and  5-7.  However  the  images  were  corrupted 
with  additive  white  noise.  The  BNR  for  both  is  10  dB.  These 
images  were,  respectively,  processed  with  the  LMS  and  RLS 
algorithms  and  the  threshold  was  set  to  11.4. 

Earlier  in  the  chapter  a  model  of  the  two  error  processes 
operating  in  a  prediction  error  image  was  presented.  It  can 
be  seen  in  Figure  5-2  that  if  the  separation  between  the  two 
processes  is  large,  the  threshold  can  be  set  to  achieve  a  high 
probability  of  target  detection  while  maintaining  a  low 
probability  of  false  alarms.   Refer  to  Tables  6  and  14  in  the 
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Figure  5-6.   Significance  Test  Results  on 
an  LMS  Processed  Image  (No  Noise) 
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Figure  5-7.  Significance  Test  Results  on 
an  RLS  Processed  Image  (No  Noise) 
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Figure   5-8.     Test  Results  on  a   LMS 
Processed  Image  (BNR=10dB) 
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Figure   5-9.     Test  Results  on  a  RLS 
Processed  Image  (BNR=10dB) 


Appendix  and  recall  the  discussion  (Chapter  III)  on  the  excess 
MSE  resulting  from  use  of  the  LMS  algorithm.  The  RLS 
algorithm  by  its  nature  does  not  suffer  from  this  type  of 
excess  MSE.  Thus  it  would  be  expected  that  the  whitened 
background  process  PDF  would  be  broader  for  the  LMS  than  for 
the  RLS.  This  result  can  be  seen  in  Table  6  and  14  by 
comparing  the   false   alarm  probabilities   at  which   100% 
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detection  occurs.  The  false  alarm  probabilities  essentially 
are  a  measure  of  the  tail  area  in  both  tails  of  the  background 
error  process.  It  is  seen  from  the  tables  that  the  broadened 
LMS  error  process  adds  an  additional  9.25%  to  the  tail  area 
(false  alarm  probability) .  In  another  set  of  results,  Tables 
3  and  11,  the  increased  false  alarm  probability  is  seen  to  be 
approximately  6.5%. 

The  effect  of  the  increased  background  error  variance  can 
also  be  seen  in  Figures  5-6  and  5-7.  In  general  from  the 
tables  it  may  be  seen  that  at  the  same  significance  level  the 
number  of  false  alarms  present  remains  approximately  constant 
for  both  LMS  and  RLS .  This  result  supports  the  proposed 
thresholding  method  and  independence  of  the  test  statistic  in 
(5-18)  .  However  the  number  of  true  targets  found  to  be  above 
the  threshold  is  smaller  for  the  LMS  method  than  for  the  RLS 
method. 

The  effect  of  the  target  error  process  in  the  impulse 
response  support  region  can  be  seen  in  Figures  5-6  and  5-7. 
The  scan  direction  for  these  figures  was  a  row-wise  raster 
pattern.  The  broadening  of  the  error  residual  in  the  vicinity 
of  the  detected  targets  can  be  seen  in  the  figures.  It  can 
be  seen  that  the  largest  error  residual  occurs  at  the 
locations  within  one  pixel  of  the  true  target  location. 

It  is  possible  to  use  the  spread  in  the  thresholds  as 
an  indicator  of  the  spread  in  the  target  process.   Bear  in 
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mind  that  comparison  between  tables  is  subject  to  the  effects 
of  squaring  the  error  in  the  test  statistic.  Also  note  in  the 
tables  the  translation  of  the  distribution  of  the  thresholds 
is  due  to  increases  in  the  target  magnitude.  These  points 
lend  support  to  the  assumed  model. 

The  preceding  comparisons  were  the  results  of  testing 
without  added  noise.  It  would  be  expected  that  adding  noise 
would  tend  to  increase  the  variance  of  the  background  error 
process  while  narrowing  the  target  error  process.  The  latter 
is  due  to  the  reduced  ability  of  the  adaptive  filter  to  follow 
the  changing  background  process.  However  in  Tables  3  through 
5,  for  the  LMS  process  with  light  noise  (BNR=30  dB)  the 
performance  is  actually  better  than  with  no  noise.  In 
adaptive  gradient  search  algorithms,  overlapping  input  data 
sets  introduces  correlation  between  the  filter  coefficient 
error  vector  and  the  data  under  the  prediction  mask.  This 
creates  an  error  bias  in  such  algorithms  as  LMS  [Ref.  13]. 
The  addition  of  a  small  amount  of  white  noise  tends  to  reduce 
the  correlation  therefore  decreasing  MSE.  This  would  decrease 
the  width  of  the  background  error  process  and  increase  the 
variance  of  the  target  process.  In  Table  5  the  beneficial 
effect  is  overcome  by  the  increased  variance  of  the  target  and 
background  processes.  Note  that  in  Tables  11  through  13  for 
the  RLS  method,  the  background  error  process  broadens  as 
expected  and  as  indicated  by  the  increased  false  alarm 
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probability.  In  Table  12  some  benefit  is  seen  in  light  noise 
(due  to  target  process  shifting) ,  but  considerably  less  than 
occurred  in  the  LMS  case.  A  visual  example  of  the  effect  of 
noise  on  the  two  methods  is  shown  in  Figures  5-8  and  5-9. 

From  the  plots  and  tables  it  can  be  seen  that  the  RLS 
algorithm  is  superior  to  the  LMS  algorithm  particularly  in  the 
area  of  dim  targets  and  in  the  presence  of  noise.  Of  course 
the  RLS  algorithm  is  computationally  more  intensive  than  the 
LMS  algorithm.  The  LMS  computational  advantage  however  may 
be  offset  by  the  additional  computational  burden  of  processing 
the  increased  number  of  false  alarms. 

A  point  of  particular  interest  is  a  region  which  may  be 
called  the  "knee"  of  the  significance  testing  process.  This 
is  the  point  at  which  the  false  alarm  rate  approximately 
doubles  for  a  10%  increase  in  detection  probability.  This 
point  appears  in  all  of  the  tables  and  occurs  at  approximately 
80%  detection  probability  for  an  LMS  whitened  image  and  90% 
detection  for  RLS  whitening.  While  these  are  not  limits  on 
the  detection  process  they  may  be  useful  as  figures  of  merit 
when  comparing  methods  of  whitening  the  image. 

D.  CONCLUSIONS 

It  has  been  proposed  that  the  decorrelation  of  the 

background  process  results  in  a  separation  of  the  target 
process  from  the  background  process.   The  experimental  data 
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supports  this  concept.  Detection  of  targets  at  low  false 
alarm  rates  requires  an  adaptive  PEF  that  produces  the 
smallest  possible  error  variance.  Thus  it  is  not  surprising 
in  the  comparison  of  the  LMS  and  RLS  methods  that  the  RLS 
results  are  better.  The  false  alarm  rates  using  this 
detection  method  are  relatively  low  for  the  detection  rates 
achieved  (see  Appendix) .  This  is  due  to  the  ability  of  the 
filter  to  adapt  to  local  changes  in  the  image.  This  in  turn 
results  in  prediction  coefficients  which  more  closely  model 
the  local  process  and  leads  to  a  more  complete  decorrelation 
of  the  image  by  the  PEF. 

The  addition  of  white  noise  to  images  degrades  the 
algorithm  performance  but  not  to  the  extent  that  might  be 
anticipated.  The  noise  broadens  the  background  error  process 
density  function  thus  adding  more  false  alarms  for  a  given 
threshold  setting.  However  the  target  error  process  variance 
is  reduced  and  offsets  some  of  the  effect.  Light  noise 
(BNR=3  0  dB)  was  found  to  increase  the  detection  rate, 
particularly  for  LMS.  Inevitably  severe  degradation  occurs 
as  the  noise  level  increases  above  10  dB. 

Significance  testing  showed  some  capability  for  detecting 
targets  with  TNR's  near  0  dB.  However  the  false  alarm  rate 
at  that  noise  level  may  be  higher  than  can  be  tolerated.  The 
conclusion  is  that  the  methods  used  show  considerable  promise 
and  capabilities  in  highlighting  dim  targets  with  TBR's  of  (1 
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to  2  dB)  and  in  the  presence  of  light  to  moderate  noise 
levels. 
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VI.  TARGET  TRACKER 

A.   GENERAL 

Throughout  this  thesis  it  has  been  assumed  that  the 
targets  to  be  located  are  pixel  or  sub  pixel  in  size.  It  has 
been  shown  that  with  the  methods  addressed  in  this  thesis  that 
target  pixels  can  be  identified.  However  this  method  also 
produces  false  alarms.  The  number  of  false  alarms  has  been 
shown  to  be  contingent  on  the  significance  level  selected. 
If  the  significance  level  is  set  to  achieve  higher 
probabilities  of  target  detection,  the  number  of  false  alarms 
increase.  Follow-on  processing  is  reguired  to  reduce  the 
number  of  false  alarms  and  to  track  the  targets  spatially  in 
successive  images. 

The  tracking  method  proposed  in  this  chapter  is  based  on 
the  multi-dimensional  pulse  tracker  proposed  by  Therrien  in 
[Ref.  26].  The  tracker  to  be  developed  links  the  two 
objectives  stated  in  the  previous  paragraph  through  the  target 
motion.  Three  possibilities  for  the  motion  of  the  suspected 
target  are  no  motion,  uncorrelated,  and  correlated.  Targets 
without  motion  will  persist  at  the  same  location  in  the  image 
through  time.  Uncorrelated  motion  is  spatial  displacement 
between  images  which  is  entirely  random  and  unpredictable. 
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Correlated  motion  is  random  spatial  displacement;  however  the 
target  motion  exhibits  a  trend  in  direction  and  rate  of 
movement . 

A  tracker  can  be  used  as  an  observer,  providing  movement 
information  that  can  be  used  to  categorize  the  targets.  This 
information  and  a  priori  knowledge  of  the  target  dynamics  can 
be  used  to  reduce  the  target  list.  This  is  accomplished  by 
comparing  the  observation  to  a  specific  profile  of  target 
dynamics.  If  the  observation  fits  the  profile,  the  target  is 
retained  in  the  target  list;  otherwise  it  is  dropped.  This 
thesis  will  not  discuss  the  logic  necessary  to  reduce  the 
target  list.  The  preceding  was  presented  to  indicate  a 
possible  means  of  using  the  tracker  as  a  additional  source  of 
target  information.  Later  in  the  chapter  it  will  be  shown 
that  some  target  dynamic  information  can  be  incorporated  into 
the  tracker  to  help  in  tracking. 

Processing  described  in  the  previous  chapters  provides  the 
location  of  possible  targets  and  access  to  two  types  of  image. 
These  images  are  the  actual  observed  image  and  the  prediction 
error  image.  In  these  images  the  target  has  a  deterministic 
shape  which  can  be  exploited  in  the  construction  of  a  tracker. 
It  is  assumed,  as  before,  that  the  target  in  the  observed 
image  is  approximately  one  pixel  in  size  and  has  a  non-zero 
magnitude.  It  was  subsequently  shown  that  the  target,  after 
whitening  with  the  PEF,  experiences  point  spreading  due  to  the 
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filter's  impulse  response.  An  image  containing  a  target  with 
intensity  profile  s(n,m)  can  be  modeled  as 

r(n,m;t)  =  s(n+f  ,m+f  ;t)  +  w(n,m) 

n      tn 

(6-1) 

where  r(n,m;t)  is  the  observed  intensity  at  a  pixel  (n,m)  in 
an  image  at  time  t.  The  term  s(n+f  ,m+C  ;t)  represents  the 
target  which  should  have  occurred  at  (n,m)  but  instead 
occurred  at  (n+f  ,m+f  )  .   The  terms  C  and  C  are  the  position 

n      m  n         m 

error  in  the  target  location  at  time  t.  This  can  be  viewed 
as  a  translation  of  the  target  from  the  expected  location  due 
to  motion.  The  term  w(n,m)  is  a  additive  noise  term.  The 
observed  image  in  the  absence  of  a  target  can  be  modeled  as 

r(n,m;t)  =  w  (n,m)  +  w(n,m) 

(6-2) 
where  w  (n,m)  is  colored  noise  (background  or  incompletely 
whitened  background)  and  w(n,m)  is  additive  white  noise.  The 
tracking  problem  is  simplified  if  the  term  w  (n,m)  is  near 
zero  or  if  the  target  power  is  considerably  larger  than  the 
colored  noise  component  power.  In  either  case  the  colored 
noise  component  can  be  ignored.  For  these  reasons  it  was 
decided  to  perform  the  tracking  in  the  prediction  error  image. 
During  the  background  whitening  stage  the  PEF  coefficients 
are  continually  adjusted  to  produce  the  minimum  error 
variance.    This  produces  a  near  minimal  phase  filter  and 
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results  in  most  of  the  impulse  response  power  occurring  near 
the  target  location.  With  these  considerations  the  target 
intensity  s(n,m)  in  (6-1)  was  modeled  as  a  spread  gaussian 


pulse  of  the  form 
r 


2 
C  exp 


s(n,m)  =  < 
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S     ^J 
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elsewhere 


v 


(6-3) 


The  term  a  controls  the  spread  of  the  modeled  signal.  For  a 
given  PEF  impulse  response  support  area,  a  is  fixed. 
Experimentally  the  best  value  was  found  to  correspond  to  be 
a  radius  around  the  target  in  which  approximately  90%  of  the 
impulse  response  power  occurs.  For  the  adaptive  3x3  PEF  this 
typically  occurred  at  a  radius  of  about  one  pixel.  The  scale 
factor  C  is  used  to  match  the  signal  magnitude  and  polarity. 
Polarity  is  positive  if  the  target  intensity  is  above  the 
local  image  mean  and  negative  if  it  is  below.  In  order  to 
match  the  model  to  the  target  signal,  an  estimate  of  the 
magnitude  and  polarity  of  the  suspected  targets  are  computed 
by  the  tracker  when  they  are  passed  by  the  detection 
algorithm. 


Although  the  true  target  intensity  profile  may  not  be 
Gaussian,  this  form  is  convenient.  It  will  be  seen  later  that 
this  tracker  is  quite  robust  with  respect  to  the  detailed 
assumptions  about  target  shape. 
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The  detection  algorithm  identifies  a  possible  target  to 
the  tracker.  The  tracker  then  extracts,  from  the  prediction 
error  image,  a  grouping  of  pixels  centered  on  the  target 
pixel.  This  will  be  referred  to  as  the  target  tracking  window 
and  its  center  point  will  be  called  the  centroid.  A  square 
tracking  window  is  used  with  the  method  developed  in  this 
thesis.  The  window  should  be  large  enough  to  include  most  of 
the  filter  impulse  response.  Since  two  targets  of  equal 
magnitude  falling  completely  in  the  tracker  window  would  be 
treated  as  if  there  were  only  one,  the  window  maximum 
dimension  is  based  on  the  minimum  distances  between  adjacent 
tracks. 

Once  the  target  has  been  identified  and  the  target 
location  has  been  passed  to  the  tracker,  track  processing 
becomes  autonomous.  The  tracker  will  compute  successive 
estimates  of  the  target  location  in  future  images.  With  each 
estimate,  the  new  tracking  window  centroid  coordinates  are 
computed  and  the  tracking  window  is  repositioned  at  the  new 
centroid.  This  autonomous  tracking  method  permits  parallel 
processing  for  detection  and  tracking  but  may  require  a 
separate  prediction  error  image  storage  buffer  for  the 
tracker.  The  buffer  is  required  if  the  processing  speed  of 
the  whitening  and  detection  processor  is  faster  than  that  of 
the  tracking  algorithm. 


79 


Additional  logic  would  be  needed  to  compare  subsequent 
targets,  identified  by  the  detection  processor,  with  existing 
tracks.  Although  not  specifically  addressed  in  this  thesis 
some  of  the  possible  features  of  the  logic  would  be:  same 
target  discrimination,  split  and  merge  of  tracks,  lost  track 
identification,  and  so  on.  Ideally  the  complete  tracking  and 
logic  package  would  produce  a  steady  state  number  of  false 
alarms  tracks.  This  occurs  when  the  number  of  new  tracks 
started  equals  the  number  of  old  tracks  terminated.  A  block 
diagram  of  the  relation  between  the  processing  blocks  in  the 
detection  and  tracking  system  is  shown  in  Figure  (6-1) . 


- 
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TU 
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Figure   6-1.     Detection  and  Tracking 
System  Block  Diagram 


The  overall  system,  less  connective  logic,  is  a  Pre-Whitening 
Unit  (PWU) ,  a  Detection  Processing  Unit  (DPU) ,  a  Tracking 
Window  Control  Unit  (WCU) ,  and  the  Tracking  Unit  (TU) .  The 
tracking  window  control  unit  and  the  tracking  unit  comprise 
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what  has  been  referred  to  as  the  "tracker."  The  WCU  controls 
the  positioning  of  the  tracking  window  based  on  the  estimates 
in  the  target  position  provided  by  the  tracking  unit.  The 
tracking  window  is  then  repositioned  so  that  the  centroid  is 
at  the  estimated  target  location  in  successive  images. 

B.   TRACKER 

In  the  derivation  of  the  tracker  to  follow  the  target 
position  error  (C    ,C   )    in  (6-1)  is  assumed  to  be  a  zero  mean 

n   in 

gaussian  random  process.  Further  it  is  assumed  that  the 
components  of  the  error  in  the  n  and  m  direction  are 
independent  and  that  the  noise  in  the  models  (6-1)  and  (6-2) 
is  a  zero-mean  white  gaussian  random  process. 

The  origin  of  the  coordinate  system  (n,m)  is  the  centroid 
of  the  tracking  window.  The  modeled  target  image  (6-1)  is 
located  at  an  assumed  position  (C  ,  C  )  which  is  to  be 
estimated.  A  diagram  of  the  tracking  window  and  its 
terminology  is  shown  in  Figure  6-2. 

Initially  an  estimator  for  the  position  error  of  a  single 
target  pixel  will  be  developed.  This  will  then  be  extended 
to  all  the  pixels  in  the  tracking  window.  To  begin  with,  some 
notation  will  be  defined.   Let  the  vector  f  be  defined  as 

'"& 

m 

(6-4) 
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Figure   6-2. 

Terminology 


Diagram   of   Tracking 


Further  let  z  be  the  vector  process  of  position  errors  in  the 
target  intensity  at  a  pixel  location  in  the  tracking  window, 
that  is 

n  n      n  i    m  m      m 

(6-5) 

Thus  the  vector  z    is  comprised  of  the  position  error  for  a 

given  pixel  in  the  tracking  window  in  each  of   k  successive 

images. 

The  vector  r  is  similarly  defined  as  the  sequence  of  image 
intensities  at  a  given  pixel  location  in  the  tracking  window, 
that  is 
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r(n,m;t*) 
r(n,m;t^) 

r(n,m;tk) 

(6-6) 
Finally  the  target  vector  s(£)  is  a  column  vector  containing 
the  sequence  of  target  intensity  values  present  at  the  same 
location  (n,m)  in  the  tracking  window,  that  is 


sfn+f^m+f^-t1) 


s(£)  = 


s(n+r^,m+rf;t2) 

n      m 


s(n+r^m+f|;;tk) 


(6-7) 
The  index  t1  ,  for  i=l,2,3...k,  is  the  time  index  for  the  k 
image  frames  considered  in  estimating  the  position  error. 
Each  pixel  intensity  is  a  sample  of  the  continuous  signal  (6- 
3)  within  the  tracking  window.  The  observed  target  is  assumed 
to  be  deterministic  but  subject  to  a  random  position  error  and 
additive  white  noise.  The  situation  is  depicted  in  Figure 
6-2. 

The  goal  of  the  tracking  algorithm  is  to  predict  the 
position  error  in  the  target  in  succeeding  images.  One  method 
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of  accomplishing  this  is  to  use  MAP  estimation.  Given  the 
statistical  models  for  the  position  error  process  and  assuming 
for  the  moment  that  the  form  of  the  target  image  process  is 
known  and  given  by  (6-3)  then  it  is  possible  to  compute  the 
a  posteriori  conditional  density  function  p(z.|r).  From  this 
density  it  is  possible  to  determine  the  best  estimate  z  in  the 
MAP  sense.  It  will  be  shown  that  future  estimates  of  £  can 
be  formed  recursively  from  this  estimate.  [Ref.  26]  Using 
Bayes  rule  the  a  posteriori  probability  density  p(z|r)  can  be 
evaluated  from 

p(zlr)  =  P(r|z)p(z) 
P(r) 

(6-8) 
Thus  using  (6-1)  the  conditional  probability  density  function 
p(r|z)    is 

p(r|z)    =  1 r79exp["  I    (r-s  (£)  )  TR  _1  (r-s  (£)  )  ] 

(6-9) 
where  R  is  the  covariance  matrix  of  the  additive  white  noise 
process.  The  a  priori  probability  density  p(z)  based  on  the 
assumptions  stated  earlier  is 

pU)  = vhr> — n72exP[  "  x  ^TV^] 

(27T)k/':|R2|1/<i  2 

(6-10) 
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where  R  is  the  covariance  matrix  of  the  position  error  in  the 
n  and  m  directions.  Because  of  the  assumed  independence  of 
£"  and  f  this  matrix  can  be  partitioned  as 


R"   ;  0 

z   I 

*2  -I— t~ 


l    J 

(6-11) 

In  order  to  find  an  MAP  estimate  we  need  to  find  the  value  of  z 

that  maximizes  (6-8) .  Since  the  denominator  is  not  a  function 

of  z    it  is  sufficient  to  find  z    to  maximize  the  numerator. 

If  we  substitute  (6-9)  and  (6-10)  into  the  numerator  of  (6- 

8) ,  take  the  logarithm,  and  drop  the  constant  term  we  are  left 

with 


C(£) 


=  -^[(r-s(£))TRw  !  (£-*<£))  +  zTRz  lz] 


(6-12) 
To  determine  z    (the  optimum  estimate  of  z)  the  derivative 

— o  p  t 

of  (6-12)  is  taken  with  respect  to  the  position  error  vector  z 
and  the  result  is  set  equal  to  the  zero  vector. 


dz 


(£)   =    ^(d    R   1  (£-£(£) 

_dz  " 


)  -  R  *Z    =    0 
'      z   — 


Thus  z    is  defined  implicitly  by 

— o  p  t 

R 


z 

— opt 


3sT(£) 


dz 


Rw  l(r-s(C)) 


(6-13) 


(6-14) 


85 


where  the  partial  derivative  term  is  a  matrix  whose  components 
are  [ds./d£.].  The  signal  present  in  any  given  image  is 
subject  to  the  position  error  in  only  that  image.  Therefore 
the  partial  derivative  matrix  has  the  form 


ds(U    = 


.1 


0 


3s  (CM 

0    BsJSh 


afm1 
0 


0 

dsJlll 
3$V 


0 
.  0 


ds(Ck) 

acnk 


.  0 


3fmk 


(6-15) 

At  this  point  it  is  beneficial  to  return  to  Equation  (6- 

14) .   Now  recall  that  the  noise,  w,  was  assumed  to  be  white. 

Therefore  the  covariance  matrix  R  is  a  diagonal  matrix  with 

components  equal  to  the  noise  power  N  .  Thus  if  we  take 
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advantage  of  this  and  the  partitioning  of  (6-15) ,  the  estimate 
for  z  can  be  written  as 


±-      (r-s(£)) 


z  (n,m)  =  RnpsT(C) 


(6-16) 


and 


z  (n,n)  =  Rm 


dsTfQ 


«f. 


^-   (r-s(£)) 

N0 


(6-17) 
where  z  (n,m)  and  z  (n,m)  are  the  position  error  estimates  in 

— n  — m  c 

the  n  and  m  directions,  respectively,  and  are  optimal  in  the 
MAP  sense. 

Now  define  vectors  a(n,m)  and  £(n,m)  as 

a(n,m)  = 


dsT(C) 


and 


£(n,m)  = 


K. 


dsT(Q 


^-   (E-s(O) 

N0 


(6-18) 


dL 


m  _J    U 


2-      (£-s(I)) 


(6-19) 
Then  the  two  estimators  can  be  written  in  the  compact  form 


z  (n,m)  =  R  a(n,m) 
_n  \  i    i  z  _ \    / 


(6-20) 


and 


zm(n,m)  =  R™  £(n,m) 


(6-21) 
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It  will  be  assumed  for  the  moment  that  the  covariance  matrices 
Rn  and  R"  are  known.  Equations  (6-20)  and  (6-21)  then  each 
represent  a  set  of  k  scalar  equations.  The  last  equations  of 
these  sets  have  the  forms 

P 


C  (n,m;p)  =  2  R"  (p-i)cr(n,m;  i) 
i=p-Tc+l 


(6-22) 


and 


rn(n,m;p)  =  2   R*  (p-j  )  |8  (n,m;i) 
i=p-k+l 


(6-23) 
where  a(n,m;i)  and  /3(n,m;i)  are  the  iin  components  of  a  and  §_ 
defined  as 


th 


a(n,m;i)  =  -L- 
N0 


dsJ(L) 


k  : 


(r-s(i)) 


(6-24) 


and 


^(n,m;i)  = 


1_  psT(£)"l 

N0  LaCJ 


(r-s(£)) 


(6-25) 
These  terms  are  the  inner  product  of  the  gradient  of  the 
signal  multiplied  by  the  difference  of  the  observed  and 
estimated  images.  These  equations  provide  the  means  by  which 
the  target  position  estimate  is  updated. 

Consider  the  form  of  s(£)  Equation  (6-7)  in  conjunction 
with  the  estimator  Equations  (6-22)  and  (6-23) .   It  can  be 
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seen  that  in  order  to  compute  the  recursive  estimates  for  the 
position  error,  k  previous  estimates  are  required.  The  upper 
limit  in  the  summation  represents  the  most  recent  image  time 
i=p=0  representing  the  beginning  of  the  track.  It  can  be  seen 
that  the  terms  of  (6-22)  and  (6-23)  with  indices  i<p  represent 
the  track  history  prior  to  time  i=p.  These  terms  for  a  causal 
system  are  identically  zero.  At  the  beginning  of  the  track 
the  observed  target  location  is  passed  to  the  tracker.  Then 
the  tracking  window  centroid  is  placed  at  the  observed  target 
location.  By  design  the  modeled  target  s(n+f  ,m+f  )  is  fixed 
at  the  center  of  the  tracking  window.  Therefore  the  estimate 
of  the  position  errors,  f  and  f  ,  at  the  beginning  of  the 
track  is  known  and  is  incorporated  in  the  centroid  position. 
Thus  all  the  essential  terms  are  available  recursively  to 
compute  the  future  estimates. 

The  terms  £"  (n,m;p)  and  f  (n,m;p)  pertain  to  a 
particular  pixel  within  the  tracking  window.  To  compute  the 
new  position  of  the  tracking  window  centroid,  these  terms  are 
averaged.  Specifically,  let  Ap  and  Ap  be  the  change  in  the 
centroid  position.  Then  the  change  in  the  centroid  position 
is  computed  from 

+l     p 

AP  +  1  =   2  Rn(p-i)a    (i) 
n  z    c  avg 

i=p-k+l 

(6-26) 
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and 


ap+1    =      2  Rm(p_i)0        (i) 

m  z  '     avg 

i=p-k+l 


(6-27) 


where  we  have  defined 

W-l  W-l 


2 


a       (i)    =   1        2  2  a(n,m;i) 

avg  wZ   n=-W+l     m=-W±l 


2  2 


(6-28) 


and 


W-l  W-l 


2  2 

0       (i)   =  1       2  2         0(n,m;i) 

avg  wZ   n=-W+l      m=-W±l 

2  2 

(6-29) 

The  results  of  (6-2  6)  and  (6-27)  represent  the  predicted  error 

in  the  target  location  in  the  next  image.   The  new  centroid 

position  is  computed  as 


p+1.   .  .  =  xp  .   .  .  +  A  p+1 

(6-30) 


x 

centroid      centroid      n 


and 


yp+1   -,  =  Yp    .,  +  A  p+1 

^  centroid     icentroi 


d 

(6-31) 

The  entire  procedure  is  repeated  in  each  succeeding  image 

frame.   In  the  absence  of  a  priori  knowledge  of  the  future 

target  position,   the  error  in  the  new  estimated  target 
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location  is  assumed  to  be  zero.  It  can  be  seen  that  the  only 
information  that  needs  to  be  retained  to  continue  the 
recursion  is  the  k  past  values  of  a   (i)  and  B        (j). 

a  v  g  a  v  g 

C.   TRACKER  DYNAMICS 

The  purpose  of  this  section  is  to  provide  some  insight 
into  the  dynamics  of  the  tracker.  Further  a  method  of 
incorporating  the  target  dynamics  into  the  design  of  the 
tracker  is  shown. 

In  the  following  analysis  the  motion  of  the  target  through 
the  sequence  of  images  will  be  assumed  to  be  uncorrelated  (see 
section  A  of  this  chapter)  .  This  type  of  motion  is  consistent 
with  the  problem  of  tracking  small  point  targets  with 
displacement  between  image  frames  of  approximately  one  to  two 
pixels.  Since  tracking  is  performed  separately  in  the 
horizontal  (n)  and  vertical  (m)  directions,  we  will  only 
consider  a  single  direction  (horizontal)  in  the  analysis.  All 
results  shown  also  apply  to  tracking  in  the  vertical 
direction. 

It  follows  from  the  assumption  of  uncorrelated  target 
motion,  that  the  error  in  the  target  position  in  the  current 
image  frame  is  independent  of  its  position  error  in  the 
previous  frame.  As  a  result,  the  position  error  correlation 
matrix  R  is  diagonal.  The  target  position  error  estimates 
for  the  p+lst  image  (6-26)  and  (6-27)  then  become 
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W-l  W-l 

2  2 

NnW^  n=-W±l  m=-W+l       ~w[ 

2  2 


'0* 


(6-32) 


and 


W-l     W-l 


i  2       2 

AP+1  =   1      2       2      3s(n+fn,m+rjr^(r-s(£)) 

NnW*   n=-W+l   m=-W+l       ^? 

2        2 

(6-33) 

where  rn(0)   and  rm(0)   are  the  values  of  the  correlation 

z  z 

function  at  lag  zero.  Since  it  was  assumed  that  the  position 
error  process  has  zero  mean,  these  terms  are  also  the  position 
error  variances.  Further  note  that  the  gradient  is  reduced 
to  a  scalar  gradient  of  the  modeled  target  image  which  is 
evaluated  at  each  point  in  the  tracking  window.  In  order  to 
simplify  the  notation  in  this  section,  the  scalar  gradient 
term  in  (6-32)  will  be  written  as  v  (n,m) . 

With  Equations  (6-1)  and  (6-32)  the  mean  of  the  position 
error  estimate  in  the  horizontal  (n)  direction  is  found  to  be 

E{A  p+1}  =  rz(0)   2  2V  (n,m)E{r(n,m;p)-s(n,m)  } 

n  o  n 

NnW^   n  m 
u     (n,m)e{W,W) 

-  rz(°)   2  S  V  (n,m)E{w(n,m;p) } 

NnW2   n  m 
u    (n,m)e{W,W) 

=  0 

(6-34) 
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where  w(n,m;p)  is  the  additive  zero  mean  white  noise  in  the 

pth  image  frame.   Further  the  variance  of  the  position  error 

estimate  is 

E{(AnP+1)2}    =   (z(0M      E{    2     S      Vn(n,m)  (r  (n,m;p) -s  (n,m;p)  ) 
\    NqVT  /  n     m 

(n,m)e{W,w} 


S   2  V  («,k)(r(*,k;p)-s(«,k;p))> 
I        k 
(2,k)e{W,W} 


Nn   /   VT  n   m 


Nq\  W  /    n   m 


(n,m)e{W,W) 
E 

NoV  W  / 

(n,m)e{W,W) 

(6-35) 

Note  that  if  the  tracker  is  to  follow  target  position 

changes   in  the   image,   the  variance   of  Ap    should  be 

proportional  to  the  variance  of  the  change  in  the  frame  to 

frame  position.   Thus  the  tracker  can  be  optimized  to  follow 

targets  with  a  particular  velocity  by  adjusting  the  term 

rn(0)  . 

z 

Figures  6-3  and  6-4  show  plots  of  the  target  position 
error  estimate  versus  the  actual  target  position  error.  In 
Figure  6-3  the  estimates  are  for  target  motion  completely  in 
the  horizontal  (n)  direction.  In  Figure  6-4  the  estimates  are 
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Figure  6-3.    Estimated  Position  Error 
Versus  Actual  Position  Error; 
a^=2.0,  r"(0)=.7  for  .5  Pixel  per 
Frame,  rn(0)=1.2  for  1.0  Pixel  per  Frame 
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Figure  6-4.  Estimated  Position  Error 
Versus  Actual  Position  Error  (Target 
Track  45  Degrees);  a  =2.0,  rn(0)=.7  for 
.5  Pixel  per  Frame,  rn(0)=.5  for  1.0 
Pixel  per  Frame 


for  target  motion  at  45  degrees  away  from  the  horizontal  with 
a  horizontal  component  of  .5  or  1.0  pixels  per  image.  In  both 
plots  it  can  be  seen  that  the  estimator  response  to  the  actual 
position  error  is  divided  into  two  regions.  The  first  is  a 
linear  region  centered  at  the  origin  and  extending  for  about 
+1  pixel  in  true  position  error.  The  remainder  is  a  region 
where  the  response  is  nonlinear. 

For  these  plots  the  term  rn(0)  has  been  adjusted  to  place 
the  known  maximum  displacement  between  frames  in  the  linear 
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region  of  the  estimator  response.  The  region  beyond  this 
point  is  the  nonlinear  region.  In  this  region  the  estimates 
of  the  target  position  error  do  not  match  the  actual  target 
position  error.  Thus,  the  result  of  setting  the  position 
error  variance  to  match  a  specific  target  motion  is  that  those 
with  greater  motion  will  tend  to  be  rejected. 

It  is  apparent  from  the  plots  that  only  a  small  portion 
of  the  tracking  window  {-3,3}  is  contained  in  the  linear 
region.  This  is  primarily  a  result  of  the  gradient  of  the 
target  model.  Error  contributions  between  the  modeled  target 
and  the  observed  target  (6-33)  receive  considerably  less 
weight  as  the  distance  from  the  centroid  increases.  A 
possible  means  of  extending  the  effective  area  of  the 
estimator  is  through  the  use  of  a  constant  magnitude  gradient. 
The  form  of  the  estimators  for  the  "constant  gradient"  method 
is 


W-l     W-l 
2        2 


V   =  l  ?    2       2      -2sgn(n)r^(r-s(0) 
NqVT  n=-W+l  m=-W+l 


(6-36) 


and 


W-l     W-l 
2        2 
Amp+  =   1  2  2       2      -2sgn(m)r™(r-s(0) 


m 


NqW"  n=-W+l  m=-W+l 


(6-37) 
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where  the  gradient  in  (6-32)  and  (6-33)  has  been  replaced  by 
-2sgn(n)  and  -2sgn(m).  Plots  of  the  estimator  response  with 
the  constant  gradient  are  shown  in  Figures  6-5  and  6-6.  The 
value  of  r  was  decreased  to  produce  estimators  with  the  same 
power  as  those  in  the  previous  plots.  The  target  motion  in 
Figures  6-5  and  6-6  is  the  same  as  that  in  Figures  6-3  and  6- 
4  respectively. 
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Figure  6-5.  Estimated  Position  Error 
Versus  Actual  Position  Error;  o  =2.0, 
rn(0)=.l  for  .5  Pixel  per  Frame,  rn^0)=.2 
for  1.0  Pixel  per  Frame 
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Figure  6-6.  Estimated  Position  Error 
Versus  Actual  Position  Error  (Target 
Track  45  Degrees);  o  =2.0,  rn(0)=.l  for 
.5  Pixel'  per  Frame,  rn=.2  for  1.0  Pixel 
per  Frame 


It  can  be  seen  that  using  a  constant  gradient  flattens  the 
response  of  the  tracker  in  the  nonlinear  region.  The  effect 
of  maintaining  the  same  power  is  to  reduce  in  the  magnitude 
of  the  estimate  near  the  origin.  However,  this  is  offset  by 
the  increased  magnitude  of  the  estimates  further  from  the 
origin.  Thus,  the  effective  response  region  has  been 
increased  virtually  to  the  window  limit  (-3,3). 
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D.   IMPLEMENTATION  OF  THE  TRACKER 

After  the  images  are  processed  by  the  target  detection 
algorithm  the  resulting  targets  and  prediction  error  image  are 
passed  to  the  tracker.  The  targets'  locations  are  placed  in 
the  tracker's  "Track  File"  until  retreived  by  the  tracking 
window  control  unit.  After  retrieval  of  the  target  data  the 
tracking  window  control  unit  accesses  the  prediction  error 
image  (see  Figure  6-7)  .  The  data  within  the  computed  7x7 
tracking  window,  centered  on  the  designated  target 
coordinates,  is  extracted. 
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Figure  6-7.   Tracker  and  Trackfile 

At  this  point  an  estimate  of  the  polarity  and  magnitude 
of  the  target  is  computed  by 


99 


3     3 
C  =  e(0,0)-  12    2   e(n,m) 
46   n=-3  m=-3 

(n,m)*{(0,0)  (0,1)  (1,0)} 

(6-38) 

where  e(n,m)   is  the  prediction  error  residue  within  the 

tracking  window.   The  term  e(0,0)  is  the  target  intensity, 

located  at  the  tracking  window  centroid.    To  compute  the 

estimate  the  intensities  of  pixels  considerably  affected  by 

the  impulse  response  due  to  the  target  have  been  excluded. 

As  mentioned  previously  these  pixels  are  generally  within  one 

pixel  of  the  target. 

After  this  computation,  both  the  tracking  window  data 
array  and  the  value  of  the  scaling  factor  C  are  passed  to  the 
tracking  unit.  The  tracking  unit  computes  an  estimate  of  the 
target  position  error  using  (6-32)  and  (6-33)  or  (6-36)  and 
(6-37) .  The  estimate  is  passed  to  the  tracking  window  control 
unit  which  computes  an  updated  tracking  window  centroid  and 
repositions  the  tracking  window.  This  new  position  and  the 
estimated  value  of  C  are  stored  in  the  track  file  until  a  new 
prediction  error  image  becomes  available  and  the  track  file 
has  been  completely  processed. 

A  block  diagram  of  the  tracking  unit  is  shown  in  Figure 
6-8.  The  summer  computes  the  difference  between  the  tracking 
window  array  and  the  target  model  array  and  the  output  is 
multiplied  by  the  reciprocal  of  the  estimate  for  the  noise 
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Figure  6-8.   Tracking  Unit  Diagram 

power.  Two  array  multipliers  are  used  to  multiply  this  result 
by  the  target  model  gradient  with  respect  to  the  n  and  m 
directions.  For  the  general  method  represented  by  (6-2  6)  and 
(6-27)  the  output  of  the  array  multiplier  is  averaged  and  the 
result  is  stored  in  a  stack.  The  averaging  is  accomplished 
by  the  block  designated  as  LPF.  The  oldest  values  of  a  (i) 
and  B       (i)  are  shifted  out  of  the  stack  and  discarded.   The 

avg 

contents  of  the  a    and  (3         buffers  are  vector  multiplied  by 

avg         avg 

the  k  lags  of  the  position  error  autocorrelation  function.  The 
result  is  the  updated  estimate  of  the  position  error.  For 
uncorrelated  target  motion  represented  by  (6-32)  and  (6-33) 
use  of  the  stack  is  not  reguired. 
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Gradient  and  modeled  target  intensity  values  for  locations 
within  the  tracking  window  are  provided  to  both  the  summer  and 
the  array  multiplier  by  the  block  labeled  "Nonlinearity  S(0 . " 
This  block  is  implemented  as  a  look-up  table. 

The  system  is  initialized  with  the  initial  target 
coordinates,  the  previously  estimated  C  term,  and  by  setting 
a       (i)=£   (i)=0  for  all  lags. 

a vg        a vg 

E.   TEST  RESULTS 

The  following  section  is  comprised  of  tracking  runs  at  a 

variety  of  noise  levels  and  for  several  different 
trajectories.  All  results  were  produced  with  the  same 
uncorrelated  trajectory  model  (i.e.,  only  the  zero-th  lag  of 
r  was  not  equal  to  zero) .  The  results  are  presented  in  Table 
6-1  and  in  Figures  6-9  through  6-20.  The  table  contains,  for 
each  figure,  the  signal  to  noise  ratio,  the  estimated  variance 
of  the  position  error  r  (0)  ,  and  the  resulting  sample  error 
variance  (Var(n)  and  Var(m))  in  the  estimates. 

In  the  first  trajectory  (Figures  6-9  through  6-14)  the 
target  motion  was  .5  pixels  per  image  frame  in  both  the  row 
and  column  direction.  In  the  second  trajectory  (Figures  6-15 
through  6-2  0)  the  target  was  moving  at  approximately  1  pixel 
per  frame.  All  figures  in  the  simulations  show  motion 
through  20  image  frames.  This  trajectory  started  at  the  lower 
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TABLE  6-1.  TRACKING  SIMULATION  RESULTS 


Figure   # 

SNR(dB) 

rz(0) 

Var(n) 

Var(m) 

6-9 (p)3 

6.4 

0.7 

.075 

.075 

6-10(p) 

2.8 

0.7 

.114 

.134 

6-ll(p) 

1.2 

0.7 

.293 

.147 

6-12(c)4 

6.4 

0.1 

.388 

.310 

6-13(c) 

2.8 

0.1 

.498 

.383 

6-14(c) 

1.2 

0.1 

.558 

.190 

6-15(p) 

6.4 

1.2 

.071 

.063 

6-16(p) 

2.8 

1.2 

.245 

.418 

6-17 (p) 

1.2 

1.2 

.323 

.277 

6-18(c) 

6.4 

0.2 

.355 

.316 

6-19(c) 

2.8 

0.2 

.513 

.698 

6-20(c) 

1.2 

0.2 

.753 

.344 

p  indicates  that  the  pulse  gradient  method  was  used, 
c  indicates  that  the  constant  gradient  method  was  used, 
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Figure  6-11.   Pulse  Gradient  Track  SNR=1.2dB 
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Figure  6-12.   Constant  Gradient  Track  SNR=6.4dB 
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Figure  6-13.   Constant  Gradient  Track  SNR=2.8dB 
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Figure  6-14.   Constant  Gradient  Track  SNR=1.2dB 
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Figure    6-16.      Pulse  Gradient  Track  SNR=2 . 8dB 
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Figure  6-17.   Pulse  Gradient  Track  SNR=1.2dB 
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Figure  6-18.   Constant  Gradient  Track  SNR=6.8dB 
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Figure  6-19.   Constant  Gradient  Track  SNR=2 . 8dB 
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Figure  6-20.   Pulse  Gradient  Track  SNR=1.2dB 
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left  corner  and  moved  to  the  upper  right  corner.  In  Figures 
6-9  through  6-14  a  180  degree  direction  reversal  was  made  by 
the  target,  at  the  location  (10,8),  followed  later  by  other 
abrupt  changes  in  the  trajectory.  The  figures  show  that  the 
estimated  trajectory  follows  the  track  well  at  all  noise 
levels.  The  table  further  shows  that  the  computed  error 
variances,  Var(n)  and  Var(m)  are  consistently  less  than  one 
pixel . 

The  estimates  of  the  position  error  were  computed  using 
(6-3) ,  (6-32)  ,  and  (6-33)  evaluated  at  the  discrete  points  in 
the  tracking  window.  This  involved  floating  point  operations 
and  round  off  of  the  final  position  error  estimate  to  the 
nearest  integer  value.  As  a  result,  quantization  noise 
becomes  a  factor  in  the  position  error  estimates.  The 
quantization  level  separation  in  the  case  of  an  image  is  one 
pixel.  Therefore  assuming  a  uniform  probability  distribution 

for  the  floating  point  result,  a  quantization  noise  variance 

2  • 
of  .083  pixel   is  expected.   It  can  be  seen  for  the  higher 

signal-to-noise  level  that  the  results  are  close  to  the 

quantization  level. 

The  noise  rejection  capability  of  the  tracking  algorithm 

using  the  pulse  gradient  form  is  seen  to  be  high.   In  Figures 

6-9  through  6-11  and  6-15  through  6-17  the  estimated  tracks 

remain  within  one  pixel  of  the  true  track.   However,  as  might 

be  expected,  the  constant  gradient  method  is  more  susceptible 
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to  noise  and  shows  deviations  from  the  track  exceeding  a  pixel 
at  even  the  higher  signal-to-noise  levels  (see  Figures  6-12 
through  6-14  and  6-18  through  6-20)  .  However  it  may  be 
required  to  use  the  constant  gradient  method  if  the  target 
motion  exceeds  2  pixels  per  image  frame.  When  the  target 
shift  between  frames  exceeds  the  width  of  the  window,  the 
track  will  be  lost.  Simulations  have  shown  that  if  only  noise 
is  contained  in  the  target  window  the  track  will  stagnate  or 
drift  very  little.  The  effective  radius  for  the  pulse 
gradient  model  was  found  to  be  approximately  1.5a  ,  where  a 
is  the  target  model  shape  factor.  This  can  be  tied  to  the 
decreased  magnitude  of  the  gradient  beyond  this  value.  For 
the  constant  gradient  method,  the  effective  radius  is 
essentially  equal  to  the  window  dimension  (W)  .  However, 
because  of  this,  the  constant  gradient  method  is  more 
susceptible  to  clutter  or  targets  crossing  the  boundaries  of 
the  tracking  window. 

F.  CONCLUSIONS 

This  chapter  was  divided  into  several  major  areas.  The 
first  was  a  general  description  of  the  interface  between  the 
detection  algorithm  and  the  tracking  algorithm.  The  other 
areas  consisted  of  a  development  of  the  tracker  equations,  a 
dynamic  analysis  of  the  tracker,  and  simulation  results. 
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The  results  of  these  and  additional  simulations  verified 
the  tracker's  expected  performance.  It  was  found  through 
simulation  that  the  gradient  of  the  proposed  target  model 
limits  the  effective  rate  at  which  targets  can  move  and  still 
be  tracked.  This  prompted  consideration  of  the  constant 
gradient  model.  In  this  model  the  modeled  target  intensity 
was  left  unchanged  but  the  gradient  was  made  a  constant  with 
appropriate  sign.  Use  of  the  constant  gradient  was  found  to 
increase  the  track  noise  and  decrease  the  resolution  (i.e., 
minimum  spacing  between  adjacent  tracks)  of  the  tracker. 

Either  form  of  the  tracker  could  be  used  as  a  target 
dynamics  observer.  An  observer  is  a  system  which  can  produce 
estimates  of  parameters  (e.g.,  velocity)  that  are  not  readily 
available.  This  observer  information  could  be  incorporated 
into  the  overall  track  file  reduction  process  proposed 
earlier.   This  may  be  a  profitable  area  for  future  research. 
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VII.  CONCLUSIONS  AND  RECOMMENDATIONS 

A.  GOALS 

The  purpose  of  this  thesis  was  to  develop  algorithms  for 
the  detection  and  tracking  of  dim  point  targets  in  infrared 
images.  Images  of  this  type  include  those  produced  by 
satellite-mounted  infrared  sensors.  The  targets  are 
restricted  to  an  approximate  size  of  one  image  pixel  and 
typically  move  at  rates  of  approximately  1  to  2  pixels  between 
successive  image  frames.  A  goal  was  also  to  develop  these 
algorithms  so  they  would  be  easy  to  implement  and  apply  to 
other  images  in  other  spectral  bands. 

B.  SOLUTION 

The  solution  presented  involves  three  steps.  The  first 
is  to  filter  the  image  with  an  adaptive  two-dimensional  fixed 
size  prediction  error  filter  (PEF) .  This  has  the  effect  of 
whitening  the  image  background  and  producing  spikes  at  points, 
such  as  targets,  different  from  the  background.  Both  two- 
dimensional  Least  Mean  Square  (LMS)  and  2-D  Recursive  Least 
Squares  (RLS)  filters  were  used.  Both  of  these  PEFs  proved 
effective  in  whitening  the  background  processes,  however,  the 
RLS  method  proved  to  be  superior  to  the  LMS  method  in  the  low 
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signal-to-noise  ratio  regime.  This  was  shown  to  be  important 
in  the  next  step  of  processing. 

The  second  step  involves  significance  testing  the  error 
residual  produced  by  inverse  filtering.  This  results  in 
suppression  of  the  error  residual  produced  by  inverse 
filtering  the  background  and  highlighting  the  targets.  The 
results  presented  showed  that  the  LMS  method  could  achieve  a 
90%  detection  probability  for  0  dB  targets  in  the  presence  of 
light  noise  (30  dB  SNR)  with  less  than  9%  false  alarms.  The 
RLS  method  produced  comparable  results  for  the  same 
parameters.  However,  with  moderate  noise  of  (10  dB  SNR) ,  a 
2  dB  target  and  a  90%  detection  probability,  the  LMS  false 
alarm  rate  increased  to  15%  while  that  of  the  RLS  increased 
only  to  10%.  The  performance  of  the  RLS  detection  algorithm 
was  significantly  improved  at  target  levels  of  approximately 
1  to  2  dB.  In  light  noise  this  resulted  in  a  90%  detection 
rate  for  an  approximately  5%  false  alarm  rate.  Both  detection 
methods  were  found  to  be  easy  to  implement  and  relatively 
inexpensive. 

The  final  step  taken  in  the  problem  solution  was  the 
development  of  a  target  tracker.  This  tracker  receives 
initial  target  location  information  from  the  previous 
detection  algorithm  and  tracks  targets  through  successive 
frames  of  data.  It  was  shown  that  a  track  can  be  maintained 
with  less  than  one  pixel  error  with  signal-to-noise  ratios 


120 


near  OdB  and  for  movement  of  1  pixel  per  image  frame. 
Additional  simulations  showed  that  track  could  also  be 
maintained  for  2  pixels  per  frame  motion  with  approximately 
one  pixel  track  error.  The  tracker  was  found  to  be  robust 
with  respect  to  the  estimates  of  the  required  parameters  and 
to  the  general  assumptions  made  in  the  development. 

C.   RECOMMENDATIONS 

The  previously  mentioned  processing  involves  images  from 
only  one  focal  plane.  Often  there  are  more  than  one  focal 
plane,  in  different  spectral  bands,  available  in  the  infrared 
sensors.  The  significance  of  a  target  in  one  band  will  differ 
from  that  in  the  other  bands.  This  will  in  turn  result  in  a 
different  PDF  for  the  target  error  process.  If  the  results 
of  significance  testing  multiple  images  in  different  bands 
can  be  merged  it  may  be  possible  to  considerably  reduce  the 
false  alarm  rate.  This  is  suggested  as  a  possible  area  for 
further  research. 

In  order  to  implement  the  proposed  method  additional  work 
needs  to  be  done  on  the  track  file  management  logic.  Some  of 
the  possible  approaches  and  considerations  were  mentioned  in 
Chapter  VI  Section  A. 
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APPENDIX 


This  appendix  contains  the  results  of  the  detection 
process  using  the  LMS  and  RLS,  algorithms.  A  detailed 
discussion  is  presented  in  Chapter  V. 


TABLE  1.  TEST  RESULTS  LMS  PROCESSED  DATA  (NO  NOISE) 


BNR(dB) 

TBR(dB) 

pf5 

Pd6 

THRESHOLD 

na 

0 

.0022 

.10 

23.33 

na 

0 

.0036 

.20 

14.16 

na 

0 

.0056 

.30 

11.96 

na 

0 

.0206 

.40 

6.64 

na 

0 

.0405 

.50 

4.94 

na 

0 

.0569 

.60 

3.96 

na 

0 

.0647 

.70 

3.63 

na 

0 

.0755 

.80 

3.28 

na 

0 

.1437 

.90 

2.14 

na 

0 

.1826 

1.00 

1.67 

Pf  is  the  ratio  of  the  sum  of  non-target  pixels  above  the 
threshold  to  the  total  number  of  pixels  tested. 


Pd  is  the  ratio  of  the  sum  of  target  pixels  above  the 
threshold  to  the  total  number  of  targets  present. 
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TABLE  2.  TEST  RESULTS  FOR  LMS  PROCESSED  DATA 


BNR(dB) 

TBR(dB) 

Pf 

Pd 

THRESHOLD 

30 

0 

.0024 

.10 

20.47 

30 

0 

.0029 

.20 

16.03 

30 

0 

.0066 

.30 

9.97 

30 

0 

.0085 

.40 

8.77 

30 

0 

.0135 

.50 

7.53 

30 

0 

.0372 

.60 

5.31 

30 

0 

.0455 

.70 

4.73 

30 

0 

.0471 

.80 

4.63 

30 

0 

.0827 

.90 

3.08 

30 

0 

.1145 

1.00 

2.43 
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TABLE  3.  TEST  RESULTS  LMS  PROCESSED  DATA 

(NO  NOISE) 


BNR(dB) 

TBR(dB) 

Pf 

Pd 

THRESHOLD 

na 

2 

.0022 

.10 

22.98 

na 

2 

.0031 

.20 

14.07 

na 

2 

.0044 

.30 

11.15 

na 

2 

.0106 

.40 

8.35 

na 

2 

.0241 

.50 

6.29 

na 

2 

.0371 

.60 

5.01 

na 

2 

.0608 

.70 

3.94 

na 

2 

.0668 

.80 

3.63 

na 

2 

.1116 

.90 

2.62 

na 

2 

.1772 

1.00 

1.72 
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TABLE  4.  TEST  RESULTS  LMS  PROCESSED  DATA 


BNR(dB) 

TBR(dB) 

Pf 

Pd 

THRESHOLD 

30 

2 

.0016 

.10 

28.36 

30 

2 

.0024 

.20 

19.20 

30 

2 

.0033 

.30 

12.99 

30 

2 

.0042 

.40 

11.33 

30 

2 

.0113 

.50 

8.10 

30 

2 

.0215 

.60 

6.56 

30 

2 

.0347 

.70 

5.65 

30 

2 

.0372 

.80 

5.44 

30 

2 

.0699 

.90 

3.78 

30 

2 

.1091 

1.00 

2.76 
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TABLE  5.   TEST  RESULTS  LMS  PROCESSED  DATA 


BNR(dB) 

TBR(dB) 

Pf 

Pd 

THRESHOLD 

10 

2 

.0052 

.10 

17.65 

10 

2 

.0079 

.20 

11.18 

10 

2 

.0125 

.30 

8.65 

10 

2 

.0213 

.40 

6.74 

10 

2 

.0266 

.50 

5.72 

10 

2 

.0412 

.60 

4.90 

10 

2 

.0545 

.70 

4.24 

10 

2 

.0892 

.80 

3.12 

10 

2 

.1557 

.90 

2.24 

10 

2 

.2491 

1.00 

1.42 
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TABLE  6.  TEST  RESULTS  LMS  PROCESSED  DATA 

(NO  NOISE) 


BNR(dB) 

TBR(dB) 

Pf 

Pd 

THRESHOLD 

na 

5 

.0014 

.10 

38.09 

na 

5 

.0020 

.20 

24.77 

na 

5 

.0024 

.30 

19.31 

na 

5 

.0031 

.40 

13.31 

na 

5 

.0040 

.50 

11.43 

na 

5 

.0132 

.60 

7.78 

na 

5 

.0269 

.70 

5.70 

na 

5 

.0472 

.80 

4.32 

na 

5 

.0621 

.90 

3.70 

na 

5 

.1529 

1.00 

2.03 
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TABLE  7.  TEST  RESULTS  LMS  PROCESSED  DATA 


BNR(dB) 

TBR(dB) 

Pf 

Pd 

THRESHOLD 

30 

5 

.0016 

.10 

38.44 

30 

5 

.0023 

.20 

22.79 

30 

5 

.0029 

.30 

16.37 

30 

5 

.0033 

.40 

13.78 

30 

5 

.0042 

.50 

11.45 

30 

5 

.0104 

.60 

8.17 

30 

5 

.0169 

.70 

7.09 

30 

5 

.0259 

.80 

6.19 

30 

5 

.0622 

.90 

4.02 

30 

5 

.1386 

1.00 

2.33 
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TABLE  8.  TEST  RESULTS  LMS  PROCESSED  DATA 


BNR(dB) 

TBR(dB) 

Pf 

Pd 

THRESHOLD 

10 

5 

.0016 

.10 

29.35 

10 

5 

.0022 

.20 

21.98 

10 

5 

.0031 

.30 

14.41 

10 

5 

.0046 

.40 

12.24 

10 

5 

.0182 

.50 

7.55 

10 

5 

.0272 

.60 

6.51 

10 

5 

.0500 

.70 

4.23 

10 

5 

.0655 

.80 

3.65 

10 

5 

.0835 

.90 

3.19 

10 

5 

.2229 

1.00 

1.79 
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TABLE  9.  TEST  RESULTS  RLS  PROCESSED  DATA  (NO  NOISE) 


BNR(dB) 

TBR(dB) 

Pf 

Pd 

THRESHOLD 

na 

0 

0 

.10 

22.99 

na 

0 

.0005 

.20 

15.06 

na 

0 

.0013 

.30 

13.21 

na 

0 

.0071 

.40 

9.74 

na 

0 

.0107 

.50 

8.42 

na 

0 

.0151 

.60 

7.59 

na 

0 

.0248 

.70 

6.70 

na 

0 

.0468 

.80 

5.20 

na 

0 

.0887 

.90 

3.44 

na 

0 

.1214 

1.00 

2.80 
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TABLE  10.  TEST  RESULTS  RLS  PROCESSED  DATA 


BNR(dB) 

TBR(dB) 

Pf 

Pd 

THRESHOLD 

30 

0 

0 

.10 

22.68 

30 

0 

.0004 

.20 

15.57 

30 

0 

.0015 

.30 

13.40 

30 

0 

.0076 

.40 

9.51 

30 

0 

.0112 

.50 

8.27 

30 

0 

.0164 

.60 

7.44 

30 

0 

.0258 

.70 

6.57 

30 

0 

.0513 

.80 

4.90 

30 

0 

.0822 

.90 

3.63 

30 

0 

.1169 

1.00 

2.88 
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TABLE  11.  TEST  RESULTS  RLS  PROCESSED  DATA 

(NO  NOISE) 


BNR(dB) 

TBR(dB) 

Pf 

Pd 

THRESHOLD 

na 

2 

0 

.10 

27.66 

na 

2 

.0001 

.20 

18.45 

na 

2 

.0003 

.30 

16.34 

na 

2 

.0035 

.40 

11.38 

na 

2 

.0054 

.50 

10.22 

na 

2 

.0066 

.60 

9.64 

na 

2 

.0133 

.70 

7.77 

na 

2 

.0234 

.80 

6.71 

na 

2 

.0579 

.90 

4.47 

na 

2 

■  ■ ...  ■ 

.1126 

—  — 

1.00 

3.03 
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TABLE  12.  TEST  RESULTS  RLS  PROCESSED  DATA 


BNR(dB) 

TBR(dB) 

Pf 

Pd 

THRESHOLD 

30 

2 

0 

.10 

27.28 

30 

2 

0 

.20 

18.47 

30 

2 

.0005 

.30 

16.25 

30 

2 

.0034 

.40 

11.51 

30 

2 

.0052 

.50 

10.00 

30 

2 

.0080 

.60 

9.27 

30 

2 

.0153 

.70 

7.39 

30 

2 

.0252 

.80 

6.46 

30 

2 

.0555 

.90 

4.69 

30 

2 

.1208 

1.00 

2.95 
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TABLE  13.  TEST  RESULTS  RLS  PROCESSED  DATA 


BNR(dB) 

TBR(dB) 

Pf 

Pd 

THRESHOLD 

10 

2 

.0005 

.10 

16.39 

10 

2 

.0051 

.20 

11.12 

10 

2 

.0120 

.30 

8.21 

10 

2 

.0255 

.40 

6.50 

10 

2 

.0427 

.50 

5.39 

10 

2 

.0519 

.60 

5.02 

10 

2 

.0711 

.70 

4.33 

10 

2 

.0917 

.80 

3.88 

10 

2 

.1085 

.90 

3.56 

10 

2 

.2157 

1.00 

2.35 
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TABLE  14.  TEST  RESULTS  RLS  PROCESSED  DATA  (NO 

NOISE) 


BNR(dB) 

TBR(dB) 

Pf 

Pd 

THRESHOLD 

na 

5 

0 

.10 

36.24 

na 

5 

0 

.20 

23.14 

na 

5 

0 

.30 

20.65 

na 

5 

.0002 

.40 

15.33 

na 

5 

.0008 

.50 

12.49 

na 

5 

.0038 

.60 

11.10 

na 

5 

.0045 

.70 

10.31 

na 

5 

.0141 

.80 

7.66 

na 

5 

.0257 

.90 

6.26 

na 

5 

.0604 

1.00 

4.49 
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TABLE  15.  TEST  RESULTS  RLS  PROCESSED  DATA 


BNR(dB) 

TBR(dB) 

Pf 

Pd 

THRESHOLD 

30 

5 

0 

.10 

34.22 

30 

5 

0 

.20 

24.77 

30 

5 

0 

.30 

22.28 

30 

5 

.0001 

.40 

17.28 

30 

5 

.0005 

.50 

14.62 

30 

5 

.0036 

.60 

13.19 

30 

5 

.0044 

.70 

11.08 

30 

5 

.0131 

.80 

8.15 

30 

5 

.0267 

.90 

6.15 

30 

5 

.0682 

1.00 

4.08 
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TABLE  16.  TEST  RESULTS  RLS  PROCESSED  DATA 


BNR(dB) 

TBR(dB) 

Pf 

Pd 

THRESHOLD 

10 

5 

0 

.10 

23.42 

10 

5 

0 

.20 

18.25 

10 

5 

.0002 

.30 

15.52 

10 

5 

.0007 

.40 

13.34 

10 

5 

.0023 

.50 

11.46 

10 

5 

.0062 

.60 

9.53 

10 

5 

.0314 

.70 

6.99 

10 

5 

.0409 

.80 

5.94 

10 

5 

.0726 

.90 

5.16 

10 

5 

.1057 

1.00 

3.69 
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